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ABSTRACT 


PL /I  procedures  have  baan  davalopad  that  uaa  finlte-differ- 
anca  techniques  to  analysa  ablation  problem*  by  embedding  than  In 
lnverse-heat-conduction  problana  with  no  moving  boundarlea.  The 
procaduraa  form  a sat  of  subroutlnas  that  can  ba  called  from  a 
problem-oriented  main  program  written  by  the  uaar.  The  procedures 
Include  provisions  for  one-,  two-,  or  three-dimensional  conduction, 
parallel  modes  of  heat  transfer,  thermal  contact,  choices  of  Im- 
plicit and  explicit  difference  techniques,  temperature-dependent 
and  directional  thermal  properties,  radiation  relief,  aerodynamic 
heating,  chemical  ablation,  and  material  removal  from  combinations 
of  flat,  cylindrical,  and  spherical  surfaces.  This  report  Is  meant 
to  serve  as  a source  of  underlying  theory  not  covered  elsewhere  and 
as  a user's  manual  for  the  PL/I  procedures.  Also  Included  are  use- 
ful debugging  elds  and  external  Identifiers,  a directory  of  Applied 
Physics  Laboratory  computer  librerles  pertelning  to  the  PL/I  pro- 
cedures, and  an  Illustrative  problem  as  an  example. 
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1.  INTRODUCTION 


The  work  described  Is  Che  result  of  en  effort  to  estebllsh 
sn  eccurete  finite-difference  (FD)  method  of  multidimenslonel  abla- 
tion analysis  that  can  be  used  economically  in  a general-purpose 
computer  program  applicable  to  reentry  ablation  problems  of  long 
duration.  Details  of  comparisons  of  various  FD  ablation  analysis 
methods  and  the  justification  for  selecting  the  one  discussed  here 
ere  given  in  Ref.  1. 

PL/1  procedures  for  incorporating  the  FD  version  of  the  em- 
bedding method  (EM)  of  ablation  analysis  (Refs.  1 end  2)  into  the 
BBS  Standard  Heat  Transfer  Program  (SHTP)  (Refs.  3 and  4)  have 
been  mitten  and  successfully  tested.  The  new  procedures  described 
here  are  collectively  called  SHTP-E  and  are  applicable  to  the  abla- 
tion of  a subliming  material  subject  to  aerodynamic  heating  during 
reentry  into  the  earth's  atmosphere. 

The  procedures  involve  two  features  that  are  new  in  SHTP-E 
and  require  some  explanation.  One  is  EM;  the  other  is  approximate 
factorization  (AF)  (Ref.  1),  an  implicit  FD  method  that  permits 
the  decomposition  of  an  N-dimensional  nodal  energy-balance  equation 
(EBE)  into  N one-dimenaional  nodal  EBE's. 


Ref.  1.  J.  D.  Randall,  "An  Investigation  of  Finite  Differ- 
ence Recession  Computation  Techniques  Applied  to  a Nonlinear  Reces- 
sion Problem,"  APL/JHU  ANSP-M-15,  Mar  1978. 

Ref.  2.  J.  D.  Randall,  "Finite  Difference  Solution  of  the 
Inverse  Heat  Conduction  Problem  and  Ablation,"  Proc.  1976  Heat 
Transfer  and  Fluid  Mechanics  Institute,  Stanford  University  Press, 
Stanford,  CA,  1976,  pp.  257-269. 

Ref.  3.  R.  K.  Frazer,  "Further  Additions  and  Improvements 
to  the  BBE  Heat  Transfer  Program,"  APL/JHU  EM-4274,  16  Jun  1969; 
Rev.  1,  29  Oct  1969;  Rev.  2,  22  Mar  1971;  Rev.  3,  29  Mar  1972; 

Rev.  4,  22  Jan  1973. 

Ref.  4.  R.  K.  Frazer,  "URLIM-A  Unified  Radome  Limitations 
Computer  Program,  Volume  2 — Users  Guide,"  APL/JHU  TG  1293B,  Apr 
1978. 
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This  report  briefly  deecrlbee  EM,  AF,  their  comblnetlon 
end  eeeocleted  terminology,  theory  not  covered  elsewhere,  the 
generel  nuaerlcel  elgorlthm  used,  end  the  use  of  eech  of  the  re- 
quired PL/ I proceduree.  Femlllerlty  of  the  ueer  with  the  SHTP 
(Refs.  3 end  4)  end  the  PL/I  (Ref.  S)  lenguege  Is  essumed.  Eech 
procedure  deecrlption  is  sterted  on  e new  pege  so  that  descriptions 
of  new  procedures  end  revisions  of  old  ones  cen  be  Inserted  from 
time  to  time. 


guage  Reference  Menual,"  IBM  File  Mo.  S360/S370-29,  Order  No. 


GC33-0009-4,  Oct  1976. 
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2.  DESCRIPTION  OF  EMBEDDING 


Figure  1 illustrates  some  of  the  basic  terminology  needed 
in  discussing  EM.  EM  is  used  to  determine  the  in-depth  thermal 
response  of  an  ablating  material,  which  Involves  moving  boundaries, 
by  embedding  the  ablation  problem  in  an  Inverse  problem  having  no 
moving  boundaries  and  no  material  removal.  The  embedding  problem 
is  properly  posed  by  appllng  a fictitious  surface-heat ing-rate 
distribution,  Qp,  to  a stationary  embedding  surface  that  covers 

the  receding  surface  associated  with  the  ablation  problem.  The 
embedding  model  used  is  of  the  embedding  problem,  not  the  ablation 
problem  that  it  contains.  The  fictitious  solid  between  the  embed- 
ding and  ablating  surfaces  is  retained  in  the  analysis.  However, 
the  physical  boundary  conditions  of  the  ablation  problem  are  used 
to  determine  the  fictitious  boundary  conditions  of  the  embedding 
problem.  As  shown  in  Section  4,  interpolation  can  be  used  to  find 
Qp  if  the  receding  surface  temperature  and  location  ("sp",  in  Fig. 

1)  are  known. 

Conductive  heat  fluxes  in  the  direction  of  ablation  (to  be 
called  the  A direction)  in  the  ablating  material  are  estimated 
using  new  temperature  (implicitly),  while  conductive  heat  fluxes 
normal  to  the  A direction  (to  be  called  the  B and  C directions) 
are  estimated  either  Implicitly  or  by  using  old  temperatures  (ex- 
plicitly). The  Implicit  A-direction  lntemodal  heat-flow  connec- 
tions form  columns  heated  on  the  embedding  surface  by  Qp,  and  each 

column  runs  normal  to  the  embedding  surface  at  least  through  the 
ablating  material.  These  columns  will  be  called  embedding  columns. 

To  begin  solving  a thermal-analysis  problem  using  EM,  one 
needs  to  know  the  distance  (e_  in  Fig.  1)  of  the  embedding  surface 

from  a fixed  reference  surface,  and  the  initial  location  of  the 
ablating  surface;  l.e.,  the  initial  value  of  sp  in  Fig.  1.  Often 

the  reference  surface,  embedding  surface,  and  initial  ablating  sur- 
face will  coincide  so  that  the  value  of  ep  will  be  zero  and  the 

initial  value  of  s_  will  be  zero.  External  heating  rates  are  known 

in  terms  of  fluxes  (Ref.  6),  but  nodal  energy  balances  are  estab- 
lished in  terms  of  total  heat-transfer  rates  (Refs.  4 and  5), 


Wave  Laser  Damage  Computer  Program,"  APL/JHU  TG  1268,  Dec  1974. 
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therefore  embedding-surface  node  areas  have  to  be  known  to  start 
the  calculations  if  the  embedding  surface  is  the  initial  ablating 
surface.  In  order  to  find  QFE,  the  current  recession  distance,  8£, 

and  the  distance,  d_  , of  each  grid  point,  1,  in  the  embedding 
Ei 

column  from  the  reference  surface  have  to  be  known  so  that  the  re- 
quired interpolation  can  be  carried  out. 

The  temperature  of  the  receding  surface  is  found  by  solving 
a surface-energy-balance  equation,  which  requires  ar  accurate  de- 
termination of  the  ablating  surface  temperature  gradient  in  the  A 
direction  and  the  local  inclination  of  the  ablating  surface  with 
respect  to  the  A direction  (Ref.  1).  To  determine  the  temperature 
gradient  properly,  the  position  of  the  ablating  surface  (s£)  with 

respect  to  the  neighboring  grid  points  (characterized  by  the  d_  ) 

1 

in  the  same  embedding  column  has  to  be  known.  For  a proper  deter- 
mination of  the  surface  inclination,  one  needs  to  know: 

1.  The  spatial  distribution  of  s£, 

2.  The  B and  C coordinates  of  the  neighboring  embedding 
columns , 

3.  The  type  of  embedding  surface  and  column  geometry 
(rectangular,  cylindrical,  or  spherical),  and 

4.  Whether  the  embedding  column  considered  lies  on  the 
edge  of  or  is  interior  to  a group  of  embedding  columns 
in  the  B or  C directions,  or  is  enclosed  in  the  B or  C 
direction  by  embedding  columns  whose  ablation  is  sym- 
metric about  it.  (Figure  2 illustrates  the  three 
situations.) 

The  required  initial  information  described  above  has  to  be 
supplied  by  the  user  and  a PL/I  procedure  has  been  written  to 
facilitate  the  process.  Once  the  information  has  been  supplied, 

EM  provides  a way  to  perform  rapid  and  accurate  ablation  computa- 
tions. The  EM  procedures  have  been  written  in  such  a way  that  the 
initial  embedding  data  need  not  be  revised  when  a calculation  is 
stopped  and  restarted.  A provision  for  reading  a second  file  of 
revised  recession  data  created  by  a previous  job  has  been  incor- 
porated into  these  procedures. 
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3.  APPROXIMATE  FACTORIZATION 


Aa  stated  in  Section  2,  lnternodal  connections  In  the  A di- 
rection, st  least  through  the  ablating  materiel,  are  implicit, 
while  lnternodal  connections  in  the  B and  C direction  may  be  either 
explicit  or  implicit.  For  any  node  that  is  treated  implicitly  in 
more  than  one  coordinate  direction,  the  method  of  AF  (Ref.  1)  is 
used  in  the  procedures  described  here  to  solve  the  nodal-energy- 
balance  equation.  AF  permits  the  decompoaltlon  of  an  N-dimansional 
Implicit  EBE  into  N one-dimensional  EBE's.  For  example,  if  a node 
has  implicit  connections  in  all  three  coordinate  dlrectlona,  there 
are  three  AF  equations:  one  in  each  of  the  directions.  A,  B,  and 
C,  respectively. 


Applied  to  conventional  heat  conduction  problems,  the  SHTP 
la  designed  to  solve  systems  of  lumped-parameter,  nodal-anergy- 
balsnce  equations,  each  of  which  has  the  following  very  general 
form: 


- if)  * •1JXTjm  - Tj+l) 

+ qJ+1J  At  - cJ(Tj+1  - Tj  ),  (1) 

where  T is  temperature,  1 denotes  the  node  of  interest,  j denotes 

1 

a neighboring  node  having  contact  with  node  i,  q^  is  the  distri- 
buted internal  heat  generation  in  node  1,  is  a total  conduc- 
tance that  can  refer  to  external  convection  or  radiation  or  inter- 
nal conduction  or  radiation,  C*  is  the  node's  thermal  capacitance 

(volumetric  specific  heat  times  node  volume),  1 denotes  time-level 
l,  t is  time,  At  ■ t1+1  - tfc,  and  0 1.  lnternodal  links 

with  a^  - 1 are  explicit,  those  with  a^  ■ 0 are  implicit.  Hybrid 

implicit-explicit  links  are  also  possible  but  will  not  be  considered 
hare.  For  brevity,  the  superscript  t will  be  dropped  and  l + 1 will 
be  replaced  by  a + in  most  of  the  following  discussion. 
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Internally,  SHTP-E  daflnaa  a vector  whose  1 — element  la 

qt  -£Vtj  - V ♦ < «> 

so  that  Eq.  1 becomes 


(3) 


The  sumawtlon  In  Eq.  2 applies  only  to  explicit  links  while  that  In 
Eq.  3 applies  only  to  Implicit  links.  Letting  conductances  for  im- 
plicit links  in  the  A,  B,  and  C directions  be  given  by  K^,  K® ^ , 

c 

and  Kjj  , Eq.  3 can  be  rearranged  to  give 


f1  + c*  (^-KlJAij  + i Kij  Aij 


pCiJAij)]Tt 


Q.At 


(4) 


where  A^  Is  a difference  operator  such  that 

vt  - < - TI  • 


(5) 


In  order  to  obtain  accurate  temperature  predictions,  At  should  be 
small,  therefore  Eq.  4 can  be  approximated  by 


U + Cl^Vl3)(1  + C1^-KljAiJ)’(1 


(6) 


Equation  6 Is  the  approximate  factorisation  of  Eq.  4.  Assuming 
consistent  node  definition,  Eq.  6 remains  consistent  with  the 
Fourier  heat-conduction  equation  If  the  A,  B,  and  C directions 
are  mutually  orthogonal. 
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By  defining  and  as 

®1-  (i+^Z<jVTt 


'i-  <1  + ^^V°1 


one  can  decompose  Eq.  6 into  the  following  three  one-dlmensionsl 
nodel  EBE's: 


- V±>  + Q1)6t  - C1(V1  - T1)1 
- Ut)At  - - Vt>, 


- T+)At  - Ci(T^  - Ut). 


Equations  7 and  8 show  that  £ and  V give  first  approximations  of 

T*  while  Eqs.  9 through  11  show  how  AF  can  be  used  in  a lumped-pa- 
rameter interpretation  applicable  to  SHTP-E.  The  AF  B and  C links 
also  form  columns,  just  as  the  A links  form  embedding  columns. 
Equation  9 is  applied  to  all  nodes  in  the  column  for  each  AF  C 
column.  The  resulting  system  of  equations  can  be  arranged  trldl- 
agonally  and  solved  rapidly  (Ref.  7)  to  give  all  the  values  in 

that  column,  because  all  the  T^  values  are  known  beforehand.  The 

process  is  repeated  for  all  the  AF  C columns  until  all  the  elements 
in  V are  found.  Similarly,  knowing  V,  U can  be  found  by  applying 
Eq.  10  to  all  of  the  nodes  in  all  of  the  AF  B columns,  after  which 
the  application  of  Eq.  11  to  all  AF  A (l.e.,  embedding)  columns 

will  give  T^\  If  there  are  no  implicit  links  in  one  or  more  of  the 
three  orthogonal  directions,  one  or  more  of  Eqs.  9 through  11  will 
become  trivial,  but  the  remaining  equations  will  still  be  proper  AF 
nodal  EBE's. 


Raf.  7.  G.  D.  Smith,  Numerical  Solution  of  Partial  Differ 
antlal  Equations.  Oxford  University  Press,  London,  1965 
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The  use  of  implicit  InCemodel  link*  hea  • beneficial  effect 
on  the  time  step,  At,  allowed  by  numerical  stability  considerations. 
Using  Eq.  2,  Eqs.  9,  10,  and  11  can  be  arranged  to  give 


Vi(Ci  + EKJjAt>  "E^jV*  - 

Vci  +EKJjATjAt  + fl^t. 


°i<°i  -Ek!juj 


At  - C^, 


+ E*JjAO  "EkJj  T^t  " CiUi* 


The  superscript  E in  Eq.  12  signifies  an  explicit  link.  A practical 
stability  criterion  can  be  obtained  from  Eqs.  12  through  14  using 
matrix  theory  and  the  method  of  row  or  column  norms  (Ref.  8).  With- 
out going  through  the  details,  the  matrix  analysis  Indicates  that 
the  sums  of  the  absolute  values  of  the  temperature  coefficients  on 
the  left-hand  sides  of  Eqs.  12  through  14,  which  contain  unknown 
quantities,  should  be  greater  than  or  equal  to  the  sums  of  the  ab- 
solute values  of  the  temperature  coefficients  on  the  right-hand 
side  of  those  equations,  which  contain  known  quantities. 

Por  brevity,  the  following  definitions  will  be  adopted  in 
performing  the  stability  analysis: 


®A,B,C  " Ci  + 2£*ij  * At* 


which  arm  the  appropriate  sums  from  the  left-hand  side  of  Eqs.  12 
through  14.  Prom  the  right-hand  side  of  Eq.  12,  there  are  two 


Stability  of  Linear 


D.  Lax  and  R.  D.  Rlchtmeyer,  "Survey  of  the 
Plnlte  Difference  Equations,"  Communications 
Vol.  9,  May  1936,  pp.  112-113. 
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possibilities  for  tho  sub,  og,  which  Is  needed : 


o 


If  C.  -UK®  At  £ 0 
1 J 13 


(16) 


and 


o 


E 


25>ijAt  - Ct>  If  C±  -^K^At  £ 0 . 


(17) 


The  ultimate  goal  In  eolvlng  Eqs.  12  through  14  la  to  find  T*. 
knowing  T;  tharafore  tha  stability  analysis  must  apply  to  tha 
collective  stability  of  solving  all  throe  aquations.  Tha  numbers 
under  tha  Inequality  signs  given  below  Indicate  tha  aquation  or 
inequality  from  which  that  particular  step  originates. 


(12) 


°E  “ 1 1 T-  • 

C 


(13) 


(18) 


°E  Ci°E 

C.  — or  1 * , 

1 °C  °C°B 


and 


(14) 


: * fborl*fi!s_ 

1 (19)  °C°B  °A° B°C 


(18) 

(19) 


(20) 


The  final  stability  criterion  Is  best  expressed  In  terms  of  sums  of 
conductances  S_  and  S_  , where 
l Ei 

[ ‘*‘^KiJ  "£K1J*  SE1  *^K1J  * (21) 

2 3 

Using  Inequality  20  and  Ignoring  0(At  ) and  0(At  ) terms,  the 
following  conservative  numerical  stability  criterion  results: 


C.  + 2S_  At  * o_. 


(22) 
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If  oE  la  glvan  by  Eq.  16,  Inequality  22  ehowe  that 

C1  + 2SI1At  * C±  • (23) 

which  la  alwaya  true  but,  according  to  Eq.  16,  la  baaad  on  having 
At  < C1/SE  , the  purely  explicit  Halt.  Equation  17  permlta  a 

higher  value  of  At  when  applied  to  Inequality  22: 

+ (Sx  - SE  )At  * 0.  (24) 


Inequality  24  la  aatlafled  unconditionally  If  Sj 

thoae  nodea  for  which  S_  > S_ 

B1 


1 


4S.  , but  for 
E1 


At  <mln 
1 


(25) 


If  Sx  >0,  the  upper  Halt  allowed  by  Inequality  25  la  larger  than 

C1^SE  * the  llait  only  on  the  explicit  llnka  and  C./(S 

1 1 Ei 

+ S_  ),  the  limit  for  no  Implicit  llnka. 

Ai 


0>  * 
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4.  COMBINING  EMBEDDING  AND  APPROXIMATE  FACTORIZATION 


EM  can  ba  combined  with  AP  by  Incorporating  tha  ambeddlng- 
aurfaca  boundary  condition  Into  Eq.  11,  giving 


[Hj'1!  - Tt>  ♦ Wn]  **  ‘ ci«i  - V' 


(26) 


whara  QpE  la  tha  flctltloua  haatlng  rata  appllad  to  embedding  node 
E (aae  Fig.  1),  and 


'El 


( O.BM 

(l,  E - 1. 


(27) 


For  each  embedding  column,  Eq.  26  la  applied  to  all  nodaa  In  tha 
column  to  form  a ayatem  of  equatlona  having  tha  aolutlon  (Refs.  1 
and  2) 


Tt  - + <>n41- 


(28) 


whara  T^g  and  T^g  are  obtained  aa  follows:  T^g  la  found  using 

U as  Initial  data  In  Eq.  11  and  no  haatlng  appllad  to  tha  embedding 
node  E,  and  It  Is  called  tha  adiabatic  thermal  response  (Raf.  1). 

T+  Is  found  by  setting  U • 0 and  applying  a heat  flow  of  1 (In 
UE1 

appropriate  units,  Btu/s,  for  example)  to  tha  embedding  node  E, 
and  It  Is  called  tha  unit  lncreattntal  thermal  response  (Raf.  1). 


Both  and  T^g  are  therefore  known  vectors,  but  at  this  point, 
T*  and  Q^g  In  Eq.  28  are  not  known.  The  following  discussion 
tells  how  to  find  Q^g  so  that  T*  can  ba  found  from  Eq.  28. 
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The  temperature,  T+  of  ttaa  ablating  aurfaca  aaaociatad 

-E 

with  embedding  aurfaca  noda  E,  la,  from  Eq.  28 


TAEa  + QFETUEa  * 


Ganarally  the  ablating  aurfaca  doaa  not  fall  on  a grid  point  and 

T*Km  and  have  to  be  found  by  interpolation.  For  conaiatency, 

accuracy  of  third  order  or  better  la  required.  In  the  PL/ I pro- 
cedure deacrlbed  hare,  third-order  Lagrange  interpolation  (Ref.  9) 
involving  three  aucceaalve  grid  point a,  I,  J,  and  K in  the  embedding 
column  near  the  ablating  aurfaca  (aee  Fig.  1)  la  uaed: 


- £ 

i-I,J,K 


where 


.t.  n 


J-l.J.K 

J^i 


a T+ 
*iXAE1, 


(*E  ~ dEV 
(dEi  " dEj‘ 


At  this  point,  Q_  and  T are  unknown  but  related  through  Eq.  29. 
FE  «E 

As  shown  below,  T is  found  from  the  ablating  surface  energy  ba- 
SE 

lance  equation  (SEBE).  Equation  29  »~111  then  give  QFE  after  which 
Eq.  28  will  give  all  the  temperatures  s.n  tha  embedding  column. 

In  general  terma,  the  SEBE  la 

^E^radin  + ^chem  + ^conv  ^redout  + ^rem  **cond  ™ 

Each  q term  is  a heat  flu*.  The  terma  qradln»  *chma>  9conv»  •»«* 

^ redout  *re  Adequately  described  by  Newman  (Ref.  6).  All  of  the 

q'a  are  nonlinear  functions  of  T . 

-E 

Raf.  9.  S.  D.  Conte.  Elementary  Numerical  Analysis.  McGraw- 
Hlll  Book  Co.,  New  York,  1965. 
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Mg  i«  the  cosine  of  the  angle  between  the  normal  to  the 

ablating  surface  and  the  A direction,  and  in  Ref.  1 it  is  shown 
that 


’ l1  + <Ve>2]  ’1/2. 


where  is  the  two-dimensional  gradient  operator  in  the  BC  plane. 

The  components  of  V^s  are  found  by  differentiating  a Lagrange 

Interpolating  polynomial  to  give  second-order  accurate  estimates. 
Referring  to  Fig.  2,  the  equation  for  3s  /3B,  where  E.  is  an  edge 

1 1 

column,  is 


- b(E1,E2)sE^  + b(E1,E2)sE^  + b^E^Sg  . 


For  3s„  /3B,  where  E-  is  an  interior  column,  the  proper  equation  is 
2 


TT  ■ b(E2"El>*E,  + b(E2"E2>V  + b<E2'E3>*E 


Finally,  for  a symmetry  column 


- 0 . 


Similar  equations  apply  for  3sE/3C.  The  b coefficients  are  given 
by 


b(Ei,EJ) 


<2»e  - £ B ) 

1 k-1,2,3  “k  , 


n <.E- 

k-1.2.3  Ej 
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In  Eq.  32,  <leond  la  the  heet  conducted  to  the  ablating  surface  by 
the  ablating  solid  and  la  given  by  (Ref.  1) 


*cond 


-cm 


(38) 


where  k(T  ) la  the  ablating  material's  thermal  conductivity  eval- 

*E 

uated  at  T*  and  the  surface  temperature  gradient  (3T/3A)+  is 
*E  ®E 

evaluated  in  the  A direction,  pointing  Into  the  ablating  material. 

A second-order  accurate  estimate  of  (3T/3A)+  is 


(fi). 


2 ciTt  • 

T T W A X 


i-I,J,K 


(39) 


where  I,  J,  and  K are  defined  in  Fig.  1 and 

ci  ” ^®E  ” 1 1 ^Ei  ~ 

jl*l  / Jl*l 


V * 


(40) 


By  combining  Eq.  28,  29,  and  39,  one  can  show  that  Eq.  38  becomes 


*cond 


/k«A  / / ^*E  _ + \ 

V -e  Ax/v"’  “• V 


(41) 


where 


AE 

U 


c T+ 
VaEI* 

U 


(42) 


I 
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”TUEs^vUE  **  positive  and  has  the  units  of  length, 

- T*  (v  /v  ) is  a temperature,  and  both  quantities  can  be  cal- 
UE8  AE  UE  ^ ^ 

culated  before  the  SEBE  is  solved.  Letting  L£  and  T£  be  given  by 


* UEs 

E ’ ' VUE 


CE  ■ TAEs  - T0E.  (^)  • r"PeCClVel>,> 


Eq.  40  becomes 


■m 


kr 


The  quantity  q is  the  energy  flux  loss  by  the  ablating  material, 
rem  . 

owing  to  the  removal  process  itself.  In  terms  of  mass  flux  m^. 

(Ref.  6),  specific  heat,  c,  and  average  temperature  change,  AT,  of 
material  ablated  during  time  At,  qre]n  is 


,cAT  . 


The  material  loss  covers  a distance  m^At/p,  where  p is  the  ablating 

material' 8 density  and  At  is  the  time  step  covered  by  the  calcula- 
tion. Therefore  near  the  ablating  surface  AT  is  given  approximately 
by 


- 25  - 


THE  JOHNS  HOSKINS  UNIVERSITY 

APPUED  PHYSICS  LABORATORY 

LAUREL  MARYLAND 


Combining  Eqs.  45  and  46  shows  that 


'rem 


(47) 


The  (pc)  has  been  inserted  because  SHTP-E  uses  (pc)  rather  than  c 
as  a basic  thermal  property. 

By  defining  the  first  five  terms  in  Eq.  32  as  q^,  that 
equation  becomes 


qmOT,(T+  ) + q .(T+  ) - 0. 
hEXT'  8g  cond  sE' 


(48) 


Equation  48  is  used  in  the  SHTP-E  procedures  if  the  embedding  and 
ablating  surfaces  do  not  coincide,  which  is  the  general  case.  Many 
problems  can  be  started  with  the  ablating  surface  lying  initially 
on  the  embedding  surface,  and  so  long  as  the  condition  persists. 


^FE  " qEXT^AE^’ 


(49) 


where  is  the  exposed  surface  area  of  embedding  surface  node  E. 

From  Eqs.  28  and  49,  the  proper  surface-temperature  equation  in 
this  case  is 


TE  “ TAEE  + qEXT^TE*TUEE  ' 


(50) 


Equation  48  or  Eq.  50,  whichever  is  appropriate,  is  solved 

by  Newton-Raphson  Iteration  for  T+  or  T*  . To  begin  the  process, 

E 

an  error  function  is  defined  as 


e(T  ) - T£  - T^g  - 

E 


(51) 


u 

u 


IJ 


u 


n 


n 

LJ 


0 


J 


0 

J 

1 
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if  th«  ablating  and  embadding  surfaces  coincide,  or 


t(T*  ) - qol  <T*  ) + ,cond  <1+  ) 


if  they  do  not.  In  either  case,  one  seeks  to  find  a number  T that 
is  a solution  of 


e(T)  - 0. 

First,  a rough  search  is  made  for  a solution  using 


_n+l  _n 

T ■ T - AT  sgn 


(*)• 


where  n and  n+1  denote  successive  iterations  at  the  new  time  level 
and 


( 1 if  x > 0 

sgn  (x)  - < 0 if  x • 0 (55) 

(-1  if  x < 0 

The  initial  value  of  AT  is  1000°R  and  an  initial  guess  for  T is 
made  from  initial  temperature  data  or  the  ablating  surface  temper- 
ature at  the  previous  time  level.  Whenever  the  solution  point  is 
passed,  AT  is  reduced  by  a factor  of  10  until 


e (Tn) 


< AT, 


after  which  Newton-Raphson  iteration, 


. Or") 


is  Invoked.  If  the  solution  starts  to  diverge,  Eq.  54  is  brought 
back  into  play  until  Inequality  56  is  satisfied,  after  which  Eq.  57 
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• used  again.  Tolaranca  limit a are  aat  on  |c(Tn)|,  Tn  itaalf, 
and  tha  maximum  nunbar  of  ltarationa.  Succaaaful  convarganca  la 

aaauaad  if  |e(Tn)|  and  |Tn  - T0"1 | fall  within  thalr  tolaranca  lim- 
lta  in  a nunbar  of  ltarationa  laaa  than,  or  aqual  to,  tha  naxlnum 
number  allowed.  If  thla  procaaa  laada  to  T < 0*R  or  T > 15  000*R, 
aolutlon  failura  ie  aaaumad  and  tha  program  axacutlon  la  atoppad. 

If  tha  ambaddlng  and  ablating  aurfacaa  coincide,  Eqa.  28  and 

49  are  uaed  to  find  all  the  T*  in  the  embedding  column  after  T+ 

haa  been  found  aa  deacrlbad  above.  If  thaaa  two  aurfacaa  do  not 
coincide,  Eqa.  28  and  29  ara  uaad  to  obtain 


r+-  T+  + (T+  - T+ 

1 “•  \w 


from  which  all  tha  T*  in  tha  embedding  column  ara  found  after  T 

1 *E 

haa  bean  found  from  Eq.  57. 

Before  computing  the  aurface  temperature  at  time  level  1+1, 

1+1 

the  A-direction  recession  distance  a_  is  estimated  from 


•r1  - •;  ♦ (£)*“ 


and  the  SEBE  la  evaluated  and  solved  for  T aa  described  above. 

*E 

During  thla  process  tha  mesa  flux  loas  m£+1  is  also  determined  the 
A~dlrection  recession  rata  than  can  be  found  from 


(£)’ 
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If  desired,  the  estimate  of  s*  can  be  refined  iteratively  by 
using  the  relation 

•r1  ■ < ♦ [.(WT.a- (*)*]*.  <«> 

where  G is  an  iteration  parameter  and  (dsg/dt)**1  is  found  from 
Eq.  60  on  the  previous  iteration.  As  shown  in  Refs.  1 and  2,  this 
process  is  fast  because  it  does  not  require  re-evaluation  of  T^g 

tue* 
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5.  PROCESSING  ABLATING-SURFACE  TEMPERATURE  AND 
PRESSURE-DEPENDENT  INFORMATION 


Aa  shown  In  Rsf.  6,  ths  Individual  ablating-surfaca  heating 
terms  in  the  SBBE  (Bq.  32)  depend  on  the  ablating-surfaca  tempera- 
ture and  the  temperature  and  pressure  of  the  gas  flowing  over  the 
ablating  surface.  The  purpose  of  Section  5 is  to  show  how  the 
ablatlng-surface  quantities,  which  are  Influenced  by  pressure  and 
ablating-surfaca  temperature,  are  treated  by  SHTP-E  procedures  in 
those  cases  where  this  treatment  differs  from  that  in  previous  SHTP 
ablation  procedures  (Ref.  6). 


The  computation  of  ablatlng-surface  characteristics  depends 
on  data  read  from  surface  thermochemical  tables  (Ref.  6).  These 
data  are  read  by  one  SHTP-E  procedure  and  used  later  by  another. 
The  nature  of  the  data  depends  on  whether  or  not  the  Lewis  number, 
Le,  of  the  gas  flowing  over  the  ablating  surface  is  unity.  In 
either  case,  the  Independent  variables  for  the  tables  read  are  gas 
pressure  (p)  and  the  ablatlon-to-diffusion  mass-loss  ratio  (S). 

The  dependent  variables  always  Include  the  ablatlng-surface  tem- 
perature and  the  gas  enthalpy  evaluated  at  this  temperature.  If 
Le  i 1,  the  gas  boundary-layer-edge  enthalpy  and  a diffusive  chemi- 
cal-heating term  are  also  read.  These  data  are  modified  so  that  « 
and  b,  where 


become  the  Independent  variables,  and  a modified  chemical  heating 
term  (Ref.  6)  is  added  to  the  tables.  Letting  x represent  any  one 
of  the  dependent  variables,  ir^  denote  a logarithmic  pressure  table 

entry,  and  b.  a logarithmic  mass-loss-ratio  table  entry. 


Previous  SHTP  ablation  procedures  (Ref.  6)  reprocess  the 
data  read  so  that  n and  ablatlng-surface  temperature  become  the  In- 
dependent variables,  since  It  is  desired  to  solve  the  SEBE  by  iter- 
ating with  the  ablatlng-surface  temperature.  During  the  reprocess- 
ing a uniform  ablatlng-surface  temperature  increment  is  used  In 
the  new  tables,  which  provides  poor  resolution  for  the  sublimation 


ablation  regime  end  unnecessarily  fine  resolution  In  the  diffusion- 
limited  ablation  regime.  Fortunately,  the  raw  surface  thermocheml- 
cal  data  are  properly  resolved,  and  the  following  discussion  de- 
scribes how  It  can  be  taken  advantage  of. 


For  each  embedding  column  at  a particular  time  level,  the 
SEBE  is  solved  for  a constant  surrounding  gas  pressure  that  Is 
known  beforehand.  By  Interpolating  on  ir,  a one-dimensional  table 
for  the  values  of  x(ir,b^)  for  all  values  of  J can  be  set  up  by 

letting 


where 


Applying  Eq.  64  to  all  of  the  dependent  variables  gives  a table 
that  coordinates  the  ablating-surface  thermochemical  data  for  the 
current  pressure.  The  choice  of  Independent  variable  in  this  table 
is  arbitrary,  but,  as  shown  above,  the  most  convenient  choice  for 
present  purposes  is  the  ablating-surface  temperature. 


The  constant-pressure-surface  thermochemical  table  just  de- 
scribed contains  data  that  are  not  valid  In  the  diffusion-kinetic 
ablation  regime  (Ref.  10).  For  Le  ■ 1,  Ref.  10  describes  a method 
of  modification  that  can  be  used  to  Insert  valid  diffusion  kinetic 
values  of  bj,  ablating-surface  gas  enthalpy,  and  the  modified  chem- 
ical-heating term  Into  the  table  for  those  tabular  ablating-surface 
temperatures  that  fall  in  the  diffusion  kinetic  regime.  This  modi- 
fication has  been  Incorporated  into  the  SHTP-E  ablation  procedure. 
To  date,  no  modification  is  available  for  Le  + 1. 


At  this  point,  one  has  a well-resolved  constant-pressure 
surface  thermochemical  table  (STT)  that  Is  useful  in  evaluating 
terms  in  the  SEBE  so  that  the  ablating-surface  temperature,  T , 

8E 

may  be  found.  Section  4 describes  how  the  evaluation  is  done  using 
Newton-Raphson  Iteration,  which  requires  the  evaluation  of  flrst- 


L.  L.  Perlni,  "Graphite  Oxidation  Modelling  In  Sub- 
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derivative  of  ell  taros  In  the  SEBE  with  respect  to  ablating-surface 
temperature.  While  previous  versions  of  the  SHTP's  ablation  pro- 
cedure properly  differentiate  nost  of  the  SEBE's  terms  with  respect 
to  ablating-surface  temperature,  those  terms  involving  the  blowing 
corrections  to  the  heat-  and  mass-transfer  Stanton  numbers  (Ref.  6) 
have  been  assumed  to  be  constant  in  the  iteration  process.  It  is 
valid  to  do  this  for  low-temperature  diffusion-limited  ablation, 
but  in  the  high-temperature  sublimation  regime,  the  blowing-corrected 
Stanton-number  variation  should  be  taken  into  account. 


The  heat-  and  mass-transfer  Stanton  numbers  are  denoted  here 
by  Cjj  and  and  their  corrections  for  blowing  are  and  C . . 

The  and  terms  depend  on  a modified  mesa-loss  ratio  6q, 

defined  aa: 


(66) 


The  funtlonal  relationships  describing  C ) end  ( SQ  ) have 

been  the  subject  of  many  experimental  and  theoretical  lnvestlga- 
ions.  Two  of  the  better  correlations  are  based  on  work  by  Mlckley 
at  el.  (Ref.  11)  (celled  the  Mlckley- Spalding  correlation)  end  Puts 
end  Bartlett  (Ref.  12). 

The  Mickley-Spalding  correlation  equations  are 


n m “ ° 

Hb  exp (2 AS  )-l 
o 

end 


(67) 


(68) 


Ref.  11.  H.  S.  Mlckley,  R.  C.  Ross,  A.  L.  Squyers,  end 
W.  E.  Stewart,  Heat,  Mesa,  end  Momentum  Transfer  for  Plow  Over  a 
Flat  Plate  with  Blowing  or  Suction.  NACA  TN  3208,  Jul  195A. 

Ref.  12.  K.  E.  Puts  and  E.  P.  Bartlett,  "Heat  Transfer  and 
Ablation  Rate  Correlations  for  Reentry  Heat  Shield  and  Hose  Tip 
Applications,"  AIAA  Preprint  No.  72-91,  Jan  1972. 
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The  X end  the  Stanton-number  ratio  (C^/Cg)  are  assumed  to  be  con- 
stants. Reference  13  recommends  X ■ 0.4  for  turbulent  flow  and 

X - 0.5  for  laminar  flow  over  the  ablating  body.  C„  and  C do  not 

n in 

depend  on  the  ablatlng-surface  temperature,  but  In  light  of  Eq.  66, 
which  contains  6,  C^,  and  C^,  they  are  functions  of  ablatlng- 

surface  temperature  and  pressure.  By  combining  Eqs.  66,  67,  and 
68,  one  can  show  that 


6. 


- fy  in  + 2*4  > • 


(69) 


Equation  69  can  be  used  to  find  d8  / dB,  and  the  STT  constructed 
from  Eq.  64  can  be  used  to  find  dB/dT  , where  T is  the  STT's  able- 

3 8 

tlng-surface-temperature  variable.  The  dC,_, /dB  and  dC  . /dB  can 

hd  o mb  o 

be  found  easily  from  Eqs.  67  and  68,  which,  after  using  the  chain 
rule,  lead  to 


^Hb 

dT 


(68) 


and 


d$ 

dT 


f- «(£)(*.)• 


(69) 


The  Puts  and  Bartlett  correlation  equations  have  the  forms 


and 


CHb  * CHF(6o> 


2XB0 

Cmb  * Si  exp(2XB  )-l 


(70) 


(71) 


Ref.  13.  ’^User's  Manual  Aerotherm  Charring  Material  Thermal 
Response  and  Ablation  Program,"  Version  3,  Aerotherm  Corporation, 
Mountain  View,  CA,  Report  UM-70-14,  Apr  1970. 
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where 


A - A($o) 


(72) 


The  functional  relationships  for  F(0q)  and  A(0q)  are  polynomials 

and  are  given  In  Ref.  6.  The  combination  of  Bqs.  66  and  71  yields 
the  relation 


0 - 


exp(2A0Q)-l 

21 


(73) 


Equation  73  gives  a simple  relationship  for  finding  0 given  8 

o 

when  In  fact,  one  desires  0Q,  given  8.  The  latter  relationship  Is 

obtainable  by  assuming  a set  of  several  values  of  Bq  and  setting 

up  a table  of  0 and  0q  values,  using  Eqs.  72  and  73.  In  the  SHTP-E 

ablation  procedure,  a table  coordinating  0,  0 , and  A Is  set  up 

o r 

before  any  ablation  calculations  are  performed.  With  a value  of 

T obtained  during  the  Iteration  process,  one  can  find  0 and 
E 

d0/dT  from  the  STT,  and  0 , A,  d0  /d0,  and  dA/d0  from  the  0-0  -A 
° O O o o 

table.  The  dC^/ dBQ  and  dC^/dO^^  can  then  be  found  from  Eqs.  70 

and  71,  and  finally  Eqs.  68  and  69  will  give  dC„. /dT  and  dC  . /dT 

Ho  s mb  s 

at  T - T+  . 

8 *E 

To  find  the  heat  loss  due  to  material  removal  (q  ) in 

# ^renr 

Eq.  32,  one  needs  to  know  the  local  mass-flux-loss,  m^,  for  the 
ablating  surface  in  the  embedding  column  with  embedding-surface 
node  E.  The  m^  has  a thermochemical  component,  , and  a mechan- 
ical erosion  component  so  that 


®E  " ®ETC  + “EMC 


(74) 
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Reference  14  give*  the  ablatlng-surface  temperature  and  eurroundlng 
gas-pressure  raagea  In  which  mechanical  erosion  la  significant  and 
shows  that  for  graphite 

Km  * Ve»V  • <”> 


For  the  Mlckley-Spaldlng  correlations.  It  Is  assumed  that 

H(T  ) - 0,  but  for  the  Puts  and  Bartlett  correlation,  the  func- 

tlonal  form  of  H(Tg  ) given  by  Ref.  14  Is  used.  mETC  Is  express- 
E 

ible  In  terms  of  the  local  surrounding  gas  boundary-layer  edge-mass 
flux  m„,  taken  parallel  to  the  ablating  surface,  as 

Vifi 


8 mGECmb 


Therefore 


“e  * 8 "GEc«b 


[l  ♦ 


dmg/dT^  at  Tg  - T 

£ 

can  be  found  directly  from  Eq.  77  using  the  STT  Eqs.  71  and  69, 
and,  when  appropriate,  the  relationship  for  H(T  ) given  by  Ref.  14. 


Ref.  14.  L.  L.  Perlnl,  "Review  of  Graphite  Ablation  Theory 
and  Experimental  Data,"  APL/JHU  ANSP-M-1,  Dec  1971. 
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6.  GENERAL  SOLUTION  PROCEDURE  SUMMARY 


Th«  general-solution  algorithm  Is  summarlzad  as  follows. 

For  sach  nods  (a.g..  node  1): 

1.  Follow  ths  consistent  nods  definition  rules  givsn  In 
Rsf.  1.  Nodal  control  surfaces  should  be  located  half- 
way between,  and  normal  to,  the  line  connecting  two  ad- 
jacent nodal  grid  points,  regardless  of  the  geometry 
used. 

2.  Specify  an  Initial  temperature. 

3.  Determine  the  nodal  capacitance  C^,  which  is  the  volu- 
metric specific  heat  (p c)^  multiplied  by  the  volume  of 
node  1. 

4.  Determine  lnteraodal  conductances  (j  » neighboring 

node),  basing  thermal  properties  on  old  temperatures 
(those  at  the  beginning  of  the  time  step).  A record  is 
kept  of  whether  pertains  to  an  Implicit  or  explicit 

connection.  If  the  connection  is  implicit,  a record  of 
the  directionality  of  the  connection  is  also  kept. 

5.  Accumulate  the  sum  of  explicit  conductances,  S_  and 

Bi 

the  sum  of  Implicit  conductances  S_  . 

6.  Determine  the  explicit  heating  rate  contribution  Q^. 

This  Includes  Internal  heat  sources,  external  heat  flows 
treated  expllcity,  and  explicit  internal  heat  flow  based 
on  explicit  lnternodal  conductances. 

For  all  the  nodes: 

1.  Find  the  minimum  stability  time  step.  At.  This  is  done 
by  comparing  S_  and  ST  for  each  node.  From  Inequall- 
Bi  *1 

tlaa  24  and  25,  If  ST  h S_  there  is  no  time  step 

h B1 


I 

I 
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or  a preset  maximum  value  (whichever  is  smaller)  Is 
used. 


2.  Determine  the  new  temperatures  of  nodes  with  no  Implicit 
connections: 


3.  Estimate  the  temperatures  of  nodes  with  B-  and  C-dlrected 
AF  connections. 


For  each  embedding  column 


1.  Estimate  the  A-direction  recession  distance,  using  the 
recession  rate  of  the  previous  time  level. 


2.  Determine  the  abla ting-surf ace  temperature 


3.  Compute  the  fictitious  heating  rate  and  find  all  the 
temperatures  in  the  embedding  column. 


4.  Find  the  new  recession  rate  In  the  A direction 


5.  (Optional)  Refine  the  estimate  of  the  recession  distance 
and  repeat  steps  2 through  4. 


mm hhh 
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7.  PL/I  PROCEDURE  DESCRIPTION 


As  usual  with  the  SHTP  (Refs.  4 and  5),  the  PL/I  procedures 
described  here  fall  Into  two  general  categories,  Initialization 
and  repetitive  thermal-analysis  computation.  The  user  writes  a 
problem-oriented  main  program  that  calls  the  procedures  described 
here,  as  they  are  needed.  The  Initialization  procedures  have  to 
precede  the  thermal -analysis  procedures.  Within  each  group,  some 
procedures  have  to  precede  others  and  the  calling  order  will  be 
Indicated  as  the  procedures  are  described. 

The  Initialization  procedures  are  designed  to  aid  the  user 
In  supplying  the  Initial  data  needed  to  perform  a satisfactory 
thermal  analysis.  The  procedures  perform  the  following  general 
tasks: 

1.  Allocating  storage  for  external  Identifiers; 

2.  Specifying  starting  and  stopping  times; 

3.  Specifying  Initial  temperature  data; 

4.  Reading  node  geometry  Information,  Including  quantitive 
geometric  characteristics  as  well  as  the  relationship 
of  each  node  to  other  nodes  with  which  It  has  direct 
links ; 

5.  Reading  embedding-surface  and  column-geometry  information; 

6.  Reading  trajectory  data; 

7.  Reading  thermal-property  data;  and 

8.  Reading  surface  thermochemical  data. 

The  procedures  for  thermal-analysis  computation  form  a loop 
and  must  follow  all  of  the  initialization  procedures.  The  proce- 
dures are  designed  to  help  the  user  perform  the  actual  repetitive 
calculations  necessary  to  determine  a thermal-response  history  of 
an  ablating  material  by  carrying  out  the  following  tasks. 

1.  Finding  node  capacitances  and  explicit  and  Implicit 
lnternodal  conductances; 
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2.  Finding  the  time  step  to  be  used; 

3.  Determining  temperatures  of  nodes  with  no  implicit 
connections; 

4.  Making  approximate  factorization  estimates  of  node 
temperatures  for  nodes  with  Implicit  B-  or  C-dlrection 
connections; 

5.  Applying  the  embedding  technique  to  each  embedding 
column; 

6.  Printing  temperature  distributions  and  ablating-surface 
Information  at  selected  time  levels;  and 

7.  Storing  data  In  data  sets  for  future  use  such  as  plot- 
ting, restarting  the  calculation  from  its  latest  time 
level,  and  saving  data  for  thermal-stress  calculations. 

The  individual  procedures  are  described  on  the  following 
pages.  Each  description  begins  with  the  appropriate  CALL  statement, 
followed  by  a list  of  the  arguments  of  each  procedure,  the  attri- 
bute of  each  argument,  and  a brief  description  of  its  purpose.  The 
argument  list  is  followed  by  a clarifying  discussion. 

INITIALIZATION  PROCEDURES 


The  initialization  procedures  follow,  each  procedure  being 
started  on  a page  by  itself  to  ensure  clarity.  These  procedures 
are  ST0RE,  SET,  READCPE,  APFAC,  READSU,  READTM,  READPR,  and  READABN. 


CALL  ST0RE  (?CAPS,  XCAPLIM,  C0DE) 


#CAPS 


FIXED  BIN.  The  number  of  nodes  In  the  FD  model 


The  number  of  nodes  plus  the  number  of 
indirect  eddresses  to  nodes  (see  dis- 
cussion of  READCPE  below) , or  e number 
et  leest  es  lerge  ss  the  lergest  indi- 
rect eddress. 


A code  lndlceting  the  type  of  analysis 
being  performed.  For  EM,  C0DE  - 5. 


ST0RE  dynamically  allocates  and  Initialises  to  zero-r.on- 
trolled  external  identifiers  used  by  several  of  the  other  proce- 
dures. ST0RE  must  precede  SET,  READCPE,  and  READSU. 
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SET 


CALL  SET  (TSTART,  TST0P,  TINIT,  0THER,  I0THER,  T0THER) ; 


TSTART 

FL0AT  BIN. 

The  time,  in  seconds,  for  thermal-anal- 
ysis computations  to  begin. 

TST0P 

FL0AT  BIN. 

The  time,  in  seconds,  for  thermal-anel- 
ysls  computations  to  end. 

TINIT 

FL0AT  BIN. 

The  initial  temperature,  in  °R,  of  all 
nodes  in  the  FD  system. 

0THER 

FIXED  BIN. 

The  number  of  nodes  whose  Initial  tem- 
peratures are  not  to  be  TINIT. 

I0THER 

(*)  FIXED  BIN.  An  erray  of  0THER  elements  specify- 
ing the  nodes  whose  initial  temperatures 
ere  not  to  be  TINIT. 

T0THER 

(*)  FL0AT  BIN.  An  array  of  0THER  elements  specify- 
ing the  temperatures,  in  *R,  of  the 

nodes  whose  Initial  temperatures  are  not 
to  be  TINIT. 

SET,  which  must  follow  ST0RE,  initialises  all  node  tempera- 
tures in  *R,  end  specifies  the  starting  and  ending  times,  in  sec- 
onds, for  the  thermal-analysis  computations.  It  is  identical  to 
the  SET  procedure  described  in  Refs.  A and  6. 

Using  T as  the  temperature  arrey,  SET  first  sets  all  tem- 
peratures equal  to  TINIT  using  T - TINIT:.  It  then  overldes  this 
specification  for  some  of  the  nodes  using  T(I0THER(I)) ; . 
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READCPE 


CALL  READCPE  (GE0FILE) ; 

GE0FILE  FILE.  The  name  of  the  file  containing  geo- 

metric data  about  the  nodes.  READCPE 
does  not  open  or  close  GE0FILE. 

READCPE  is  similar  to  READCP  (Refs.  4 and  6),  but  contains 
features  relating  to  EM  and  AF.  Its  purpose  is  to  read  and  store 
the  contents  of  GE0FILE,  which  contains  data  characterizing  each 
node 'a  geometry  and  its  relationship  to  other  nodes. 

The  first  card  in  GE0FILE  should  contain  a descriptive  title 
of  the  file's  contents;  the  second  card  should  have  title  headings 
for  the  remaining  cards,  which  require  the  following  information. 


COLUMNS  DATA  DESCRIPTION 

1-4  Fixed  format  (no  decimals).  A node  number  (direct 

address)  or  an  Indirect  address  (greater  than  the 
largest  node  number).  For  further  discussion,  see 
"Direct  and  Indirect  Addressing." 

3 

10-19  E format.  The  node  volume  (ft  ) for  directly  ad- 

dressed Information  or  0.E0  for  indirectly  addressed 
information. 

20-29  F format  (decimal  required).  Material  code.  This 

can  vary  from  1.  to  999999.  For  further  discussion 
see  "Material  Code  Conventions." 

30-39  F format  (decimal  required).  Number  of  neighboring 

nodes  with  directly  addressed  connections,  if  equal 
to  0.,  1.,  or  2.;  reference  to  Indirect  address  if 
greater  than  two. 

40-49  F format  (decimal  required).  First  neighboring  node. 

If  positive,  this  signifies  a conductive  link  or  a 
contact  link  with  another  material.  If  negative,  it 
signifies  a radiation  link.  If  an  A,  B,  or  C appears 
in  this  field  to  the  left  of  the  numerical  entry,  a 
type  A,  B,  or  C implicit  link  is  signified.  If  only 
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50-59 


60-69, 

70-79 


the  numerical  entry  appears,  the  link  la  explicit. 

For  more  discussion  of  the  use  of  the  A,  B,  and  C 
codes,  see  "Direct  and  Indirect  Addressing." 

E format.  The  geometric  part  of  the  conductance  to 
the  first  neighboring  node.  For  a conductive  link, 
this  field  gives  the  internodal  contact  area  in 
aquare  feet  divided  by  the  distance  In  feet  between 
the  grid  points  of  the  node  of  Interest  and  Its  first 
neighbor.  For  a contact  or  radiation  link,  It  gives 
the  appropriate  area  in  square  feet. 

These  entries  have  the  seme  significance  for  the 
second  neighboring  node  that  the  entries  in  columns 
40-59  have  for  the  firat  neighboring  node.  They 
are  used  If  the  entry  in  columns  20-29  Is  greater 
than  99,  or  if  the  entry  in  columns  30-39  Is  greater 
than  or  equal  to  two.  Zero  or  blank  entries  are 
permissible. 


The  last  card  in  GE0FILE  should  be  a blank  card.  READCPE 
must  precede  APFAC  and  follow  ST0RE.  The  numerical  entries  In  the 
seven  fields  in  columns  10-79  are  stored  in  an  external  Identifier 
called  XCAP.  For  node  I,  XCAP  (1,1)  Is  the  node  volume,  etc. 
XCAP's  declaration  Is 

DCL  XCAP  (*,7)  FL0AT  BIN  CTL  EXT; 

Users  may  find  this  Information  useful  when  defining  distributed 
heat-source  terms  that  require  node  volume. 

Direct  and  Indirect  Addressing 

If  N is  the  number  of  nodes  in  the  FD  model,  the  entry  in  columns 
1-4  is  I,  and  I * N,  then  I is  a node  number  and  the  XCAP  informa- 
tion on  the  card  pertains  directly  to  node  I.  This  is  called  di- 
rect addressing.  If  the  entry  in  columns  30-39  is  M and  M > N, 
then  further  XCAP  information  is  expected  on  a card  containing  M 
in  columns  1-4.  XCAP  (M,l)  should  be  zero,  but  the  remaining  en- 
tries, XCAP  (M,2)  - XCAP  (M,7)  will  contain  more  information  about 
node  I.  This  is  called  indirect  addressing.  If  necessary,  indi- 
rect addresses  can  be  used  to  access  more  information  on  other  in- 
directly addressed  cards.  Indirect  addressing  is  convenient  for 
connecting  more  than  two  neighbors  to  e node  and  for  specifying 
parallel  heat-transfer  links,  such  as  gas  conduction  and  radiation 
across  a gap. 
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Type  A (i.e.,  the  A direction  In  Fig.  1)  implicit  lnteraodel 
connections  must  be  initiated  by  direct  addressing  but  indirect  ad- 
dressing can  be  used  to  construct  parallel  connections.  Type  B or 
C implicit  connections  can  be  initiated  by  either  direct  or  indi- 
rect addressing.  Regardless  of  direct  or  indirect  addressing, 
implicit  links  must  always  have  an  A,  B,  or  C code,  even  if  there 
are  parallel  connections.  For  each  node  that  has  implicit  links, 

A,  B,  and  C neighbors  can  be  addressed  simultaneously,  but  not  more 
than  one  neighbor  of  each  type  can  be  addressed  directly  or  indi- 
rectly. Parallel  implicit  connections  have  to  be  made  on  separate 
cards. 

As  la  the  case  with  READCP  (Refs.  4 and  6),  nodal  connections 
should  be  specified  in  only  one  direction  in  GEdFILE.  If  node  (9 
has  been  designated  as  node  l's  neighbor  in  the  table,  the  latter 
is  understood  when  node  69  is  identified  as  node  l's  neighbor. 

Materiel  Cod*  Conventions 

Table  1 describes  the  material  code  conventions  used  in  the 
SHTP-E  procedures.  1 is  the  node  of  Interest,  N1  and  M2  are  its 
directly  or  indirectly  addressed  first  and  second  neighbors  in  col- 
umn 40-49  and  60-69,  and  a,  b,  c,  d,  e,  and  f are  any  of  the  sin- 
gle digits  (0-9).  The  table  describes  the  general  convention  and 
gives  some  examples. 

As  indicated  in  the  discussion  under  READPR,  each  individual 
material  has  a single  material  code  that  can  vary  from  1 to  99. 

Gaps  and  contacts,  as  well  as  materials  with  finite  volume,  can  be 
treated  as  individual  materials.  For  matarlals  with  finite  volume 
and  directional  conductivities,  separate  codes  can  be  defined  for 
each  direction.  In  contrast  to  READCP  (Refs.  4 and  6),  the  limit 
on  the  mother  of  materials  is  99  instead  of  9 and  no  special  mate- 
rial code  convention  is  adopted  for  contacts.  The  potentially 
harmful  effect  of  high  contact  conductance  on  the  allowable  time 
step  can  be  eliminated  by  making  the  contact  link  implicit. 


- 44  - 


I 

THC  JOHNS  HOPKINS  UNIVCNSITV 

IWED  PHYSICS  LABORATORY 

lAUMl  MAHYCANO 


Table  1 

Material  Code  Conventions  Used  in  SHTP-E 


I 
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Amc 


CALL  APPAC  (C0DE,  AFFILE); 

C0DE  CHAR(*).  C0DE  must  be  'B'  or  'C'.  It  is  used  to 

indicate  whether  the  contents  of  APFILE 
refer  to  iapllclt  lnternodal  links  in  the 
B or  C direction. 

APFILE  FILE.  The  contents  of  APFILE  contain  information 

about  the  starting  nodes  in  type  B or  C AF 
columns.  APFAC  opens  and  closes  APFILE. 

APFAC  defines  each  AF  column  (see  discussion  following  Eq. 

11  in  Section  3)  in  a way  that  can  be  understood  by  the  appropriate 
procedures  in  the  thermal-analysis  computation  loop.  READCPE  only 
reads  information  that  may  result  in  implicit  B or  C AF  lnternodal 
links;  APFAC  finalises  those  links. 

The  first  card  in  APFILE  is  a title  card;  the  remaining 
cards  contain  a list  of  the  node  numbers  of  the  starting  nodes  in 
the  AF  columns  desired.  These  entries  are  LIST  controlled,  l.e., 
format  free,  and  the  numbers  read  should  be  first  nodes  in  the  AF 
columns  read  by  READCPE.  If  APFAC  does  not  define  an  AF  column 
that  has  been  read  by  READCPE,  the  corresponding  lnternodal  links 
are  assumed  to  be  explicit,  even  though  they  have  been  tagged  as 
implicit  by  a B or  C in  READCPE.  This  is  accomplished  by  not  en- 
tering the  starting  node  number  in  APFILE.  If  APFILE  is  undefined 
(no  DD  card),  empty,  or  contains  no  starting  node  numbers,  there 
are  no  AF  links  in  the  C0DE  direction. 

There  are  also  various  checks  built  into  APFAC.  Program 
execution  will  be  terminated  if  C0DE  is  not  'B'  or  'C*  or  if  an 
attempt  is  made  to  cross  an  AF  column  with  itself  or  another  AF 
column  of  the  same  type.  However,  cyclic  AF  columns,  in  which  the 
last  node  in  the  column  is  implicitly  connected  to  the  first  node, 
are  permitted.  These  could  occur  in  the  angular  direction  around 
a sphere  or  cylinder.  APFAC  also  checks  the  AF  column  structure 
read  by  READCPE  to  make  sure  that  no  more  than  ona  implicit  B or  C 
neighbor  is  addressed  directly  or  indirectly  by  each  node  that  has 
AF  connections.  The  AF  columns  generated  by  APFAC  are  printed  out 
so  that  the  user  can  check  them  to  be  sure  they  are  constructed  as 
desired. 

APFAC  must  follow  READCPE. 
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CALL  READSU  (NB0DY,  NNPERS,  SRFDATA) 


The  total  number  of  embedding  surface 
body  stations.  A body  station  is  a 
collection  of  embedding  surface  nodes 
whose  associated  ablating  surfaces  are 
exposed  to  the  same  external  environment 


NB0DY 


(*)  FIXED  BIH.  An  array  of  NB0DY  numbers  giving 
the  number  of  embedding  surface  nodes 
In  each  body  station.  NNPERS (I)  - J 
means  that  there  are  J embedding  sur- 
face nodes  In  the  I—  body  station. 


NNPERS 


The  file  containing  embedding  and  initial 
receding  surface  geometry  Information. 
This  file  should  not  be  called  SURFACE. 
READSU  does  not  open  or  close  SRFDATA. 


SRFDATA 


READSU  is  used  to  read  embedding-surface  data,  column- 
geometrical  data,  and  related  Initial  ablatlng-aurface  data.  The 
file  SRFDATA  is  composed  of  cards  that  have  to  be  arranged  as  de- 
scribed below. 


Cards  1-6:  Text  describing  Information  on  the  following 
cards,  e.g.,  comments,  titles,  and  column  headings. 


Card  7:  An  asterisk  In  column  1 


Card  8:  The  name  of  a file  (starting  In  column  1)  contain- 
ing the  distances  in  Inches  of  the  grid  points  In  the  embedding 
column  In  the  ablating  material  from  the  fixed  reference  surface 
In  the  A direction,  i.e.,  the  d's  in  Fig.  1.  This  file  should  con- 
tain nothing  but  the  d's  In  ascending  order.  The  input  data  are 
format  free  (LIST  controlled).  The  d's  given  should  not  go  beyond 
the  back  face  of  the  ablating  matarial,  even  if  the  embedding  col- 
umn  does. 


Card  9:  A code  indicating  the  type  of  geometry  through 
which  the  embedding  columns  pass.  This  code  should  start  In  column 
1.  The  following  four  codes  are  acceptable: 


r 
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RECT:  Rectangular,  or  Cartesian,  geometry.  The  embedding 
columns  can  run  In  any  one  of  the  three  coordinate  directions. 

CYLZ:  Cylindrical  geometry.  Embedding  column  runs  In  the 
axial  direction. 

CYLR:  Cylindrical  geometry.  Embedding  column  runs  In  the 
negative  radial  direction. 

SPH:  Spherical  geometry.  Embedding  column  runs  In  the 
negative  radial  direction. 

The  next  several  cards  contain  Information  about  embedding 
surface  nodes.  The  first  two  columns  of  each  card  starting  a new 
data  group  contain  Information  about  the  nature  of  the  embedding 
surface  node  as  viewed  in  the  B and  C directions,  the  first  per- 
taining to  the  B direction  and  the  second  to  the  C direction.  Each 
data  group  describes  a single  embedding  surface  node.  An  S in  col- 
umn 1 indicates  that  the  embedding  surface  node's  column  is  a sym- 
metry column  in  the  B direction,  as  it  is  defined  in  Fig.  2.  An 
E In  column  1 means  that  the  embedding  surface  node's  column  Is  an 
edge  column  when  viewed  in  the  B direction.  A blank  or  any  charac- 
ter other  than  S or  E in  column  1 means  that  the  embedding  surface 
node's  column  is  an  Interior  column  when  viewed  in  the  B direction. 
Similar  remarks  apply  about  C direction  codes  in  column  2.  The  re- 
maining eleven  input  items  are  LIST  directed  and  may  appear  on  the 
first  card  in  the  data  group  or  following  cards  in  the  same  data 
group.  These  items  have  the  roles  described  below  and  must  appear 
in  the  order  Indicated. 

1.  Embedding  surface  node  number.  If  this  number  is  posi- 
tive, printing  of  the  temperatures  of  the  nodes  in  the  fictitious 
ablating  material  (see  Fig.  1)  is  suppressed;  if  it  is  negative, 
temperatures  for  these  nodes  are  printed. 

2.  Body  station  number. 

3.  Node  number  within  the  body  station. 

2 

4.  Embedding  surface  area  (ft  ). 

5.  A coordinate  of  the  embedding  surface  node. 

6.  B coordinate  of  the  embedding  surface  node. 

7.  C coordinate  of  the  embedding  surface  node. 


J 
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11.  The  distance,  in  inches,  of  the  embedding  surface  node's 
grid  point  from  the  reference  surface  in  the  A direction, 
i.e.,  e in  Fig.  1. 


If  embedding  surface  node  69  is  the  third  node  in  the  eighth 
body  station,  the  first  three  items  described  above  should  be 
69  8 3. 

If  there  are  no  B-  or  no  C-direction  items,  zeroes  should 
be  entered  In  the  appropriate  places. 

Each  embedding  distance  file,  giving  the  d's  in  Fig.  1,  re- 
quires a DD  card  in  the  job's  control  language.  When  a new  embed- 
ding distance  file  or  embedding  geometry  type  is  required,  the  next 
card  should  contain  only  an  asterisk,  repeating  card  7,  and  the 
above  sequence  of  data  entries  following  card  7 should  be  repeated. 

The  significance  of  the  A,  B,  and  C coordinates  mentioned 
above  depends  on  the  type  of  embedding  geometry  under  consideration, 
as  shown  in  Table  2.  The  notation  used  in  this  table  is  defined 
in  Fig.  3. 


Table  2 

Definitions  of  A,  B,  and  C 


Embedding 

Geometry 

A 

B 

C 

RECT 

x (in.) 

y (in.) 

z (in.) 

CYLZ 

z (in.) 

r (in.) 

6 (deg.) 

CYLR 

r (in.) 

z (in.) 

0 (deg.) 

SPH 

R (in.) 

♦ (deg.) 

0 (deg.) 
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The  data  suppllad  In  tha  flla  SRFDATA  pertain  to  a situation 
bafora  any  ablation  has  takan  placs,  when  all  recession  rates  are 
sero.  It  is  not  necessary  to  revise  SRFDATA  to  restart  a calcula- 
tion after  some  ablation  has  occurred  end  is  proceeding  at  a finite 
rate.  The  appropriate  data  revisions  sre  performed  by  ABAER0E,  de- 
scribed later. 

READSU  should  follow  ST0RE. 
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CALL  RZADTM  (I,  IND,  Dl,  D2,  D3,  D4,  D5,  NR,  INFILE); 


1 

FIXED  BIN. 

An  index  that  should  be  different  each 
time  READTM  is  called. 

IND 

(*)  FL0AT  BIN.  An  array  of  Independent  variables 
with  N + 1 elements. 

Dl, 

D4, 

D2,  D3, 

D5 

(*)  FL0AT  BIN.  Arrays  of  dependent  variables  each 
having  N elements. 

NR 

FIXED  BIN. 

The  number  of  ltsms  reed  from  the  arrays 
described  above,  usually  0 or  1. 

INFILE 

FILE. 

The  file  con tain! eg  the  list  to  be  read. 

READTM  (Refs.  4 and  5)  is  a general  cable  tc*4ar.  IND(l)  ■ N, 
the  number  of  table  values.  D1(J),  D2(J),  etc.  eorree*3*d  Co  DID 
(J+l). 


When  used  In  conjunction  with  SHTP-E,  READTM' e primary  pur- 
pose is  to  store  the  addresses  of  time-dependant  trajectory  data 
for  each  body  station.  I represents  the  body  station;  DID  is  a 
list  of  flight  times  (seconds);  and  Dl,  D2,  D3,  D4,  and  D5  are 
lists  of  recovery  enthalpy  (Btu/lbm),  Incident  radiant  heat  flux 
2 

(Btu/s-ft  ),  enthalpy-based  cold-wall  heat  transfer  coefficient 
2 

(lbo/s-ft  ),  local  total  pressure  (ataospherea) , and  a dummy 
variable,  respectively.  These  data  are  read  or  calculated  from 
raw  data  produced  by  the  three-degree-of-freedom  trajectory  code 
(Ref.  15)  prior  to  calling  READTM.  Therefore,  NR  ■ 0 and  INFILE  - 
SYS IN;  l.e.,  READTM  does  not  read  the  data  but  merely  stores  the 
addresses  of  the  data. 

READTM  can  also  be  used  in  conjunction  with  other  tables, 
but  care  should  be  taken  that  I is  not  one  of  the  body-station 
numbers . If  data  are  being  read  by  READTM,  the  first  two  cards  in 
INFILE  are  comment  cards  and  the  remaining  NR  cards  are  arranged 
as  described  below,  where  0 < NR. 

Ref.  15.  L.  L.  Per ini,  "User*a  Manual  for  the  3D0F  Trajectory 
Computer  Program,”  APL/JHU  ANSP-M-6,  Sep  1973. 
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DATA  DESCRIPTION 

IND  (J+l) 

DK(J) , K“l,  2,  3,  4,  or  5 (F-format,  decimals 

required) 

READTM  can  be  used  anywhere  In  Che  main  program  ahead  of  the 
thermal-analysis  computation  loop. 


COLUMNS 

1-9 

10K-10K+9 
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READPR 


CALL  READPR  (MATFILE,  ABMAT,  GAPM) ; 

MATFILE  FILE.  A file  containing  the  names  of  the  ther- 

mal-property files.  READPR  opens  and 
doses  this  file. 

ABMAT  FIXED  BIN.  The  material  code  for  the  ablating 

material. 

GAPM  (*)  FIXED  BIN.  A list  of  gap  material  codes. 

READPR  is  used  to  read  thermal  property  tables.  One  call 
to  READPR  replaces  several  calls  to  READRK  (Refs.  4 and  6). 

The  first  card  in  MATFILE  is  a comment  card.  The  remaining 
cards  contain  a LIST-directed  list  of  thermal-property  file  names 
(quotes  are  required  around  each  name).  The  first  file  listed  per- 
tains to  material  1,  the  second  to  material  2,  and  so  on.  A DD 
card  is  required  in  the  job's  JCL  for  each  file  listed.  As  many 
as  99  file  names  are  allowed. 

READPR  counts  the  number  of  material  files  in  MATFILE  and 
then  allocates  storage  to  the  thermal-property  pointer  to  cover 
that  many  materials.  It  then  opens  each  property  file  and  counts 
the  number  of  cards  in  the  file,  after  which  it  closes  the  file, 
allocates  storage  for  the  thermal  properties,  reopens  the  property 
file,  and  calls  READRKD,  a version  of  READRK  (Refs.  4 and  6)  with 
dynamically  allocated  property  pointers.  READRKD  reads  the  pro- 
perty data  and  stores  their  addresses.  When  control  is  returned 
to  READPR.  the  property  file  is  closed. 

Each  thermal-property  file  begins  with  two  comment  cards, 
followed  by  several  cards  giving  property  data  as  described  below. 

COLUMNS  DATA  DESCRIPTION 

(All  F format,  decimals  required,  blank  fields  are 
zeroes) 

1-14  Temperature  (°R) 

15-29  Volumetric  specific  heat  (Btu/ft^-aR) 

30-44  Thermal  conductivity  (Btu/h-ft-°F)  or 

o 

Contact  conductance  (Btu/h-ft  -°F) 
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45-59 

60-69 


Emisaivity , for  use  in  computing  gap  radiation 
Exposed  surface  absorptivity. 


6APM  is  a list  of  material  property  codes  that  are  intended 
to  refer  to  gaps.  The  property  files  for  these  materials  are  ex- 
pected to  list  gas  conductivities  and  radiant  Interchange  factors 
(in  the  emissivity  column)  in  the  gaps.  If  GAPM  (I)  ■ +J  > 0,  then 
the  property  data  for  material  J is  unaltered  after  it  is  read.  If 
GAPM  (I)  ■ -J,  then  the  thermal  conductivity  data  items  for  material 
J are  all  set  equal  to  zero  and  the  gap  medium  is  treated  as  a 
vacuum. 


READPR  should  precede  READABN  (described  next) 


~7.  " 
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READABN 


CALL  READABN  (#PL,  #BETA,  LEW#,  RH0,  C0NV,  AIRPRNT,  THRM0CH) ; 


#PL 

FIXED  BIN. 

The  number  of  pressure  entries  in  the 
thermochemical  tables  (Ref.  6). 

#BETA 

FIXED  BIN. 

The  number  of  ablation-to-dlf fusion  mass- 
loss  ratios  in  the  thermochemical  tables. 

LEW# 

FL0AT  BIN. 

The  Lewis  number  of  the  fluid  surround- 
ing the  ablating  surface,  usually  assumed 
to  be  1. 

RH0 

FL0AT  BIN. 

The  density  of  the  ablating  material  in 
lbm/ft1 * 3. 

C0NV 

FL0AT  BIN. 

A conversion  factor  for  converting  the 
units  in  thermochemical  table  entries 
to  English  units.  The  tables  are  usu- 
ally given  in  terms  of  calories,  grams, 
and  °K,  so  that  C0NV  is  usually  1.8. 

AIRPRNT 

FIXED  BIN. 

If  AIRPRNT  - 1,  READABN  prints  a listing 
of  the  surface  thermochemical  tables. 

If  AIRPRNT  # 1,  the  surface  thermochem- 
ical  tables  are  not  printed. 

THRM0CH 

FILE. 

The  name  of  the  file  containing  the  sur- 
face thermochemical  tables.  READABN 
opens  and  closes  THRM0CH. 

READABN  is  used  to  read  the  surface  thermochemical  tables, 
which  are  described  in  detail  by  Newman  (Ref.  6)  and  to  construct 
a table  of  ablating  material  enthalpy  as  a function  of  temperature, 
also  described  by  Newman.  READABN  is  similar  to  Newman's  READAB, 
but  has  five  fewer  arguments  because: 

1.  The  list  of  the  ablating  material's  temperstures  snd 

volumetric  specific  hests  needed  by  READABN  is  passed 

to  it  internally  by  READPR. 
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2.  Uniform-temperature- Increment  tables  are  not  uaad  In 
EM' a ablation  procadura,  ABAER0E.  Tablaa  interpolated 
by  pressure  are  constructed  from  the  original  nonuniform 
temperature  increment  tables  by  ABAER0E,  as  described  in 
Section  5. 

3.  AIRPRNT  is  now  en  argument  rather  than  an  external  Iden- 
tifier. 


READ ABN  must  follow  READER 
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THERMAL-ANALYSIS  COMPUTATION  LOOP 

The  remaining  procedures  perform  the  repetitive  tharmal- 
enelyeie  computations,  beginning  vlth  SET's  time  TSTART  end  ending 
with  TST0P.  During  the  first  pass  through  this  loop,  no  computa- 
tions are  made,  but  several  Initializations  and  storage  allocations 
are  carried  out  and  the  final  structuring  of  the  embedding  columns 
is  performed. 

The  first  procedure  called  is  CAPC0NE.  It  may  be  desirable 
to  precede  CAPC0NE  with  a definition  of  distributed  internal-heat- 
source  terms.  For  example,  suppose  that  nodes  I1  through  I2  have 

a heat  generation  rate  of  HINT  Btu/h-ft3.  To  Include  the  effect 
of  HINT  in  his  model,  the  user  should  first  write  the  declare 
statement 

DCL(Q(*) , XCAP(*,7))  FLOAT  BIN  CTL  EXT; 

in  his  main  program.  Then  the  computation  loop  would  have  the 
following  general  appearance. 

S0URCE : D0  I - II  T0  12; 

Q (I)  - XCAP(I.l)  * HINT; 

END  S0URCE; 

CALL  CAPC0NE; 


other  call  statements 


G0  T0  S0URCE; 

This  looks  like  an  infinite  loop,  but  the  called  procedures  keep  a 
record  of  the  current  time  variable  and  CAPC0NE  stops  the  execution 
of  the  program  when  the  time  passes  TST0P. 

Unless  otherwise  indicated,  the  thermal  analysis  computation 
loop  (TACL)  procedures  must  be  called  in  the  order  in  which  they 
are  described  on  the  following  pages. 


mm  ■■ 
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CALL  CAPC0KE; 

CAPC0NE  la  one  of  the  most  Important  of  all  the  SHTP-E  pro- 
ceduraa.  CAPC0NE  la  a veralon  of  CAPC0N  (Refs.  4 and  5)  adapted 
for  Implicit  AF  lnternodal  connections  and  the  new  material  code 
conventions  discussed  under  READCPE. 

On  the  first  pass  through  the  TACL,  CAPC0NE  finalizes  the 
embedding  column  structure  read  by  READCPE.  The  embedding  column 
nodes  are  printed  out  so  that  the  user  may  check  to  see  if  the  de- 
sired columns  have  Indeed  been  formed.  On  subsequent  TACL  passes, 
CAPC0NE  computes  the  thermal  capacitances  of  each  node,  lnternodal 
conductances,  sums  of  Implicit  and  explicit  conductances  for  each 
node,  and  explicit  heat  flows  between  neighboring  nodes,  all  of 
which  It  does  while  It  keeps  track  of  the  proper  material  codes 
and  the  directionality  of  implicit  lnternodal  connections.  When 
TIMS,  the  SHTP-E  time  variable,  passes  TST0P , CAPC0NE  stops  program 
execution. 

Unlike  CAPC0N  (Refs.  4 and  6),  CAPC0NE  does  not  call  a pro- 
cedure called  C0MC0N  to  compute  the  effects  of  contact  conductances. 
Contact  conductances  are  treated  as  described  under  READCPE  and 
READPR. 


CAPC0NE  must  precede  STEPE. 
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U 

STEPE 


i 

CALL  STEPE  (TSWITCH,  MINSTP,  MAX S TP) ; 


TSWITCH 

(*) 

FL0AT  BIN. 

A list  of  #TSW  times  (seconds)  at 
which  the  minimum  and  maximum  allow- 

U 

able  time  steps  should  be  changed. 

n 

MINSTP 

(*) 

FL0AT  BIN. 

A list  of  #TSW  minimum  allowable 
time  steps  (seconds) . 

MAXSTP 

(*) 

FL0AT  BIN. 

A list  of  #TSW  maximum  allowable 
time  steps  (seconds) . 

[j 

During  each  pass  (except  the  first  one)  In  the  TACL,  STEPE 
coeiputes  the  stability  time  step,  advances  TIMS  by  one  step,  com- 
putes the  temperatures  of  nodes  with  no  implicit  internodal  links, 
and  calls  a procedure  called  BAPSUB,  which  makes  B and  C direction 
AF  temperature  estimations  for  those  nodes  that  have  B and  C AF 
connections. 

The  three  arguments  passed  to  STEPE  should  have  the  same 
array  length,  say  #TSW.  TSWITCH  should  be  given  in  ascending  order; 
i.e.,  TSWITCH (1)  < TSWITCH(2)  < TSWITCH(3),  etc.  The  minimum  and 
maximum  allowable  time  steps  used  by  STEPE  are  selected  from  the 
MINSTP  and  MAXSTP  lists,  depending  on  the  current  value  of  TIMS, 
whose  units  are  also  seconds.  The  table  below  shows  how  the  se- 
lection is  made. 


TIMS 

IsL. 

< - TSWITCH(l) 

> TSWITCH (I  - 1) 
& < - TSWITCH (I) 

(I  > 1) 

> TSWITCH (#TSW) 


Minimum 
Time  Step 

Lsl 

MINSTP (1) 
MINSTP (I) 

MINSTP (#TSW) 


Maximum 
Time  Step 

, fad 

MAXSTP (1) 
MAXSTP (I) 

MAXSTP  (iPTSW) 


ij 


1 


.] 

.1 

3 

i 

1 
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If  there  Is  no  limit  on  the  stability  time  step,  or  If  the 
stability  time  step  Is  greater  than  the  maximum  allowable  time 
step,  or  If  the  maximum  allowable  time  step  Is  given  as  a negative 
number,  the  time  step  used  Is  the  absolute  value  of  the  maximum 
time  step.  If  the  stability  time  step  Is  less  than  the  minimum 
allowable  time  step,  the  program  Is  stopped.  Under  all  other  cir- 
cumstances, the  time  step  used  is  the  stability  time  step. 

STEPS  must  follow  CAPC0NE  and  precede  ABAER0E. 
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abaKe 


CALL  ABAER0E  (ID,  CMH,  LEW#,  VIEW,  TSPACE,  LAMBDA,  RH0,  TLIM,  QLIM, 
ITLIM,  PRNTM,  PLTC0DE,  ADIABATIC,  SITLIM,  G); 


ID 

FIXED  BIN. 

Body  station  number. 

FL0AT  BIN. 

Stanton  number  ratio  (C^/C^  in  Section  5) 

LEW# 

FL0AT  BIN. 

Lewis  number. 

VIEW 

FL0AT  BIN. 

View  factor  for  radiation  incident  on 
the  ableting  surface. 

TSPACE 

FL0AT  BIN. 

Space  temperature  (*R). 

LAMBDA 

FL0AT  BIN. 

Blowing  correction  factor. 

RH0 

FL0AT  BIN. 

3 

Ablating  material  density  (lbm/ft  ). 

TLIM 

FL0AT  BIN. 

Convergence  tolerance  limit  on  surface 
temperatures  (*R). 

QLIM 

FL0AT  BIN. 

Convergence  tolerance  limit  on  surface 

2 

heat  flux  imbalance  (Btu/ft  -s). 

ITLIM 

FIXED  BIN. 

Maximum  number  of  iterations  allowed  in 
finding  the  ablating  surface  temperature. 

PRNTM 

(*)  FL0AT  BIN.  List  of  times  (seconds)  at  which  de- 
tailed information  about  ablation  calcu- 
lations is  to  be  printed. 

PLTC0DE 

FIXED  BIN. 

Plotter  code.  If  PLTC0DE  # 1,  nothing 

is  saved  for  plotting.  If  PLTC0DE  - 1, 
data  ara  savad  for  the  Rublnsteln-Kleln 
plottar  coda  (Ref.  16). 


Ref.  16.  N.  Rubinstein  and  L.  E.  Klein,  "GENPLOT-A  Uear 
Oriented  Executable  Program  for  Linear,  Semi-Log  and  Log-Log  CALCOMP 
Plote  - User 'a  Guide,"  APL/JHU  FlC(2)-75-U-035,  19  Dec  1975. 
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ADIABATIC 

SITLIM 


G 


FIXED  BIN.  If  ADIABATIC  - 1,  the  receding  surface 
Is  Insulated,  otherwise  It  Is  not. 


FIXED  BIN.  The  number  of  iterations  desired  to  Im- 
prove the  estimation  of  the  current  re- 
cession distance.  For  no  Iterations, 
SITLIM  - 0. 

FL0AT  BIN.  Ablating  surface  location  Iteration 

parameter.  See  Eq.  61  and  the  follow- 
ing discussion. 


ABAER0E  solves  the  ablating  surface  energy  balance  equation 
(SEBE)  for  the  ablating  surface  temperature,  and  calculates  the 
local  recession  rate  and  distance  in  the  A direction  for  sublima- 
tion ablation.  It  also  calculates  temperatures  for  the  nodes  In 
each  embedding  column  for  both  explicit  and  internodal  links  in  the 
AF  B and  C direction.  ABAER0E  has  to  be  called  for  each  body 
station.  It  analyzes  all  the  embedding  columns'  ablating  surfaces 
for  each  body  station.  Embedding  surface  node  numbers  and  areas 
are  obtained  internally  from  data  stored  by  READSU. 

ABAER0E  computes  ablation  using  the  Mickley-Spalding  corre- 
lation (Refs.  1 and  6)  if  CMH  > 0 or  the  Putz  and  Bartlett  correla- 
tion (Refs.  6 and  13)  plus  a mechanical  erosion  modification  (Refs. 

6 and  14),  if  CMH  £ 0. 

In  many  applications,  the  use  of  a Lewis  number  of  1 is 
adequate,  but  one  can  also  use  nonunity  Lewis  numbers  in  ABAER0E. 
Care  should  be  taken  that  the  proper  data  have  been  read  by  READABN. 

The  incident-radiation  view  factor  is  usually  1,  and  for 
reentry  work  the  apace  temperature  is  generally  assumed  to  be  0°R. 

The  numerical  value  of  LAMBDA,  the  blowing  correction  factor, 
passed  to  ABAER0E  is  used  only  if  CMH  > 0.  Recommended  values  are 
0.5  for  laminar  flow  over  the  ablating  surface  and  0.4  for  turbulent 
flow.  When  CMH  ^ 0,  ABAER0E  computes  LAMBDA  internally.  See  Sec- 
tion 5 for  further  discussion. 

The  ablating  material's  density  is  assumed  to  be  constant, 
and  all  accounting  for  thermochemical  phenomena  is  confined  to  the 
ablating  surface. 

The  SEBE  solving  process  in  ABAER0E  tries  to  force  the  net 
heat  flux  imbalance  on  the  ablating  surface  to  be  zero,  taking  QLIM 
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as  an  acceptably  small  Imbalance  and  maintaining  a small  error  In 
the  ablating-surface  temperature,  lass  than  or  equal  to  TLIM,  with- 
in ITLIH  Iterations. 

Numerical  values  In  the  PBNTM  should  be  in  ascending  order 
(i.e.,  PRNTM(l)  < PRNTM(2)  < PRNTM(3),  etc.).  If  PRNTM(l)  < 0, 
ABAER0E  prints  Information  at  every  time  level.  If  all  values  of 
PRNTH  are  positive,  ABAER0E  prints  Information  only  for  each  of 
them  until  TIMS,  the  current  time,  la  greater  than  PRNTM'a  maximum 
value,  after  which  ABAER0E  prints  information  for  all  subaequent 
time  levels. 

If  PLTC0DE  ■ 1,  ABAER0E  stores  current  values  of  19  ablating- 
surface  characteristics  for  each  embedding  column  In  a file  called 
ABPL0T.  The  file  ABPL0T  should  be  allocated  to  a sequential  card 
data  set.  The  Rubinsteln-Klein  plotter  program  (Ref.  16)  can  be 
used  to  plot  the  data  saved  In  AEPL0T ’ s data  set  as  a function  of 
time.  The  19  data  Items,  in  the  order  of  their  listing  in  the  file, 
are: 

1.  Current  ablating-surface  temperature  (*R); 

2 

2.  Total  mass  loss  per  unit  area  (lbm/ft  ) In  the  embedding 
column; 

2 

3.  Current  mass  loss  rate  (lbm/s-ft  ); 

4.  Integrated  ablating  aurface  energy  Imbalance  (Btu/ft2); 

3.  Current  ablating  surface  energy  Imbalance  (Btu/s-ft2); 

6.  Current  recession  rate  (in./s)  in  the  A direction; 

7.  Total  recession  distance  (in.)  In  the  A direction; 

2 

8.  Convective  heating  rate  (Btu/s-ft  ) (for  further  expla- 
nation see  T(IPL0T+7),  page  55  of  Rof.  6); 

2 

9.  Integrated  convective  heating  (Btu/ft  ); 

10.  Heating  due  to  ablating  surface  chemical  reaction 

2 

(Btu/s-ft  ) (for  further  explanation  see  T(IPL0T+9), 
page  56  of  Ref.  6); 

2 

11.  Integrated  chemical  heating  (Btu/ft  ); 

12.  Radiative  heat  flux  absorbed  by  ablating  surface 
(Btu/s-ft2); 
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13.  Integrated  radiative  heating  absorbed  by  ablating  sur- 
face (Btu/ft2); 

14.  Radiation  relief  to  space  from  ablating  surface 
(Btu/s-ft2); 

15.  Integrated  radiation  relief  to  space  from  ablating  sur- 
face (Btu/ft2); 

16.  Conductive  heat  flux  from  ablating  material  to  ablating 

2 

surface  (Btu/s-ft  ) ; 

17.  Integrated  conductive  heat  flux  from  ablating  material 

2 

to  ablating  surface  (Btu/ft  ); 

18.  Net  external  heating  loss  given  up  to  material  removal 

2 

(Btu/s-ft  );  and 

2 

19.  Integrated  material  removal  heat  loss  (Btu/ft  ). 

The  times  at  which  these  data  are  to  be  plotted  are  read  from  a 
file  called  ABPLTM,  a card  file  containing  a title  card  followed 
by  cards  containing  a LIST-dlrected  list  of  times  (seconds).  After 
TIMS  passes  the  last  time  in  this  file,  no  further  Information  is 
passed  to  the  file  ABPL0T . 

ADIABATIC  - 1 Imposes  an  adiabatic  boundary  condition  on 
the  ablating  surface.  If  this  happens  when  the  embedding  and 
ablating  surfaces  coincide,  q£XT  is  set  equal  to  zero  in  Eq.  50. 

+ * 

If  the  two  surfaces  do  not  coincide,  T is  set  equal  to  T_  in 

SF  . L 

Eq.  44. 

Using  Eq.  61,  the  estimate  of  the  position  of  the  ablating 
surface  In  the  embedding  column  can  be  refined  by  going  through 
SITLIM  Iterations  using  G as  an  Iteration  parameter.  Equation  59 
Is  used  to  obtain  the  first  estimate  of  the  ablating  surface  posi- 
tion. If  SITLIM  " 0,  Eq.  61  Is  not  used  st  all.  For  fast  ablation 
G ■ 0.5  gives  good  results,  while  G ■ 1.5  works  better  for  slow 
ablation  (Refs.  1 and  2). 

On  the  first  pass  through  the  TACL,  ABAER0E  allocates  stor- 
age for  several  Identifiers  and,  if  CMH  * - 0,  sets  up  the  8,  Bq, 

and  X(6q)  table  discussed  In  Section  5.  Subeequently,  for  each 
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body  station  at  a given  time  level,  all  ablating-surface  energy- 
balance  calculations  are  performed  at  a constant  surrounding  gas 
pressure,  which  ABAER0E  finds  from  the  trajectory  tables  whose  con- 
tents can  be  found  from  the  addresses  saved  by  READTM.  For  the 
fixed  pressure,  ABAER0E  sets  up  the  constant-pressure  ablating  sur- 
face thermochemical  tables  described  in  Section  5.  For  Le  ■ 1, 
ABAER0E  calls  DIKI*,  a procedure  that  corrects  the  table  in  the 
diffusion  kinetic  ablation  regime.  There  is  no  correction  proce- 
dure for  nonunity  Lewis  number. 

DIKI  performs  the  required  table  corrections  using  several 
default  parameters  that  can  be  overridden  using  data  Introduced 
in  a file  called  L0WABF.  If  there  is  no  DD  card  for  L0WABF,  the 
default  table  corrections  are  made.  If  L0WABF  has  a DD  card  and 
the  file  contains  only  the  entry 

ITB  ■ 0; 

the  default  table  corrections  are  made.  On  the  other  hand  if 
ITB  • (any  negative  number); 

appears  in  L0WABF,  no  table  corrections  are  t^ade.  For  values  of 
ITB  that  are  neither  negative  nor  sero,  the  following  additional 
data  are  expected  by  DIKI  in  the  L0WABF  file.  The  first  group  of 
data  is  DATA  directed. 

3/2  1/2 

K0:  Reaction  rate  constant  (lbm/ft  s atm  ). 

E:  Activation  energy  (cal/mole). 

ZD:  Diffusion  limited  CO/CO^  ratio. 

TL0W:  Lowest  temperature  (*R)  at  which  table  corrections  are 
to  be  made. 

THIGH:  Highest  temperature  (°R)  at  which  table  corrections 
are  to  be  made. 

The  semicolon  should  follow  the  THIGH  assignment.  Following  these 
data  is  a LIST-directed  list  of  ITB  pairs  of  wall  temperatures  (°R) 
and  CO/ CO ^ ratios  in  ascending  order  of  wall  temperature  (Table  3). 

The  first  wall  temperature  should  be  TL0W  and  the  ITB^  one  should 
be  THIGH.  Reference  10  shows  how  to  select  the  parameters  to  over- 
ride DIKI  defaults,  listed  in  Table  4. 


* Private  communication  with  L.  L.  Per ini,  APL/JHU,  1977. 


j 


I 

I 

I 
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Table  3 

Wall  Temperature  and  CO/CO2  Ratio 


Wall 

Temperature 

<°R> 

CO/CO, 

Ratio 

1800 

0.001 

1980 

0.245 

2070 

0.945 

2160 

1.910 

2250 

4.560 

2430 

12.520 

2610 

16.580 

2880 

117.000 

3600 

200.000 

Table  4 

Default  Parameters  for  DIKI 
(Graphite  Values) 

K0  * 9.65  x 10  5 lbm/ft^^s  atm^^ 
E « 4.4  x 10^  cal/mole 
ZD  = 17 
TL0W  - 1800°R 
THIGH  3600°R 
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ABAER0E  sets  up  s thermochemical  table  as  described  above 

I and  treats  wall  or  gas  enthalpy,  mass-loss  ratio,  and  chemical 

heating  as  functions  of  wall  temperature. 

In  the  proceas  of  preparing  to  solve  and  solving  the  SEBE, 
ABAER0E  has  to  search  three  different  types  of  tables:  those  hav- 
ing surrounding  gas  pressure,  ablating  surface  temperature,  and 
(for  CMH  *■  0)  ablatlon-to-dlffuslon  mass-loss  ratio  as  independent 
variables.  In  any  of  these  cases,  the  Independent  variable  may 
fall  outside  the  range  of  table  values.  If  this  happens,  a message 
is  printed  but  the  calculations  continue.  If  the  Independent  vari- 
able is  less  than  the  lowest  table  value,  the  values  of  the  depen- 
dent variable  corresponding  to  the  lowest  Independent  variable  are 
used.  On  the  other  hand,  for  independent  variables  higher  than  the 
highest  table  value,  a linear  extrapolation  from  the  end  of  the 
table  is  used  to  find  the  dependent  variable. 

Once  the  constant-pressure  surface-thermochemical  table  is 
set  up,  ABAER0E  solves  the  SEBE  for  the  ablatlng-surface  tempera- 
ture for  each  embedding  column  passing  through  the  body  station. 

In  each  embedding  column,  ABAER0E  calls  BACKSUB,  a procedure  that 
computes  the  adiabatic  and  unit  incremental  thermal-response  vec- 
tors for  that  column.  After  this  is  done,  ABAER0E  calls  SURFACE, 
a procedure  that  computes  the  local  surface  inclination  with  re- 

* # 

spect  to  the  A direction,  L_  and  T_  for  use  in  Eq.  44,  and  the 

c b 

adiabatic  and  unit  Incremental  temperatures  on  the  ablating  surface. 
These  tasks  are  performed  using  geometry  information  obtained  from 
READSU.  Now  the  Newton-Raphson  solution  of  either  Eq.  51  or  52  is 
carried  out  to  find  the  ablatlng-surface  temperature.  Once  it  is 
found,  ABAER0E  computes  the  fictitious  heating  rate  on  the  embed- 
ding surface  node  and  uses  Eq.  29  to  find  the  temperatures  of  all 
the  nodes  in  the  embedding  column. 


If  one  desires  to  stop  the  job  or  Job  step  using  the  SHTP-E 
procedures  and  to  restart  the  ablation  calculation  on  a subsequent 
job  or  Job  step,  ABAER0E  csn  save  ablatlng-surface  heat-balance 
data  needed  to  do  this  at  the  end  of  one  Job  or  job  step  and  read 
it  at  the  beginning  of  a subsequent  job  or  Job  step.  Output  data 
are  saved  in  a file  called  RECRST0.  If  there  is  no  DD  card  for 
this  file,  the  file  is  ignored.  RECRST0  should  be  allocated  to  a 
sequential  data  set.  Input  data  are  read  from  a file  called 
RECRSTI,  consisting  of  data  saved  by  a previous  job's  or  job  step's 
RECRST0  file.  As  with  RECRSTI,  if  there  is  no  DD  card  for  RECRST0, 
the  file  is  ignored.  The  contents  of  either  file  are  cards  or  card 
images  (logical  record  length  of  80).  The  first  card  in  the  file 
gives  the  time  at  which  the  Job  producing  the  cards  stopped,  and 
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the  remaining  cards  contain  ablatlng-surface  heat-balance  data  in 
groups  of  four  cards  for  each  embedding  column.  This  information 
is  useful  for  setting  aside  space  for  a data  set  on  the  RECRST0  DD 
card. 

File  Summary 

There  are  five  files  that  can  be  used  with  ABAER0E,  although 
it  is  not  necessary  to  use  any  of  them.  If  PLTC0DE  ■ 1,  the  user 
has  to  supply  DD  cards  for  the  files  ABPLTM  (input)  and  ABPL0T 
(output).  The  default  diffusion  kinetic  modification  to  the  abla- 
ting surface  thermochemical  tables  can  be  overridden  using  the  file 
L0WABF.  If  there  is  no  DD  card  for  L0WABF,  DIKI  is  executed  using 
default  graphite  parameters.  Data  can  be  saved  for  restarting  the 
thermal  analysis  simply  by  supplying  a DD  card  for  the  file  RECRST0, 
an  output  file.  If  there  is  no  DD  card  for  this  file,  no  attempt 
is  made  to  save  any  ablating  surface  data  for  future  use.  Data 
from  a previous  job's  or  Job  step's  RECRST0  file  can  be  read  by 
ABAER0E  through  the  file  RECRSTI . No  attempt  is  made  to  read  such 
data  if  there  is  no  DD  card  for  RECRSTI.  The  output  files  ABPL0T 
and  RECRST0  should  be  allocated  to  sequential,  not  partitioned, 
data  sets.  All  of  the  files  are  card  files;  i.e.,  their  logical 
record  length  is  80. 

ABAER0E  must  follow  STEPE  and  precede  WRITEE  and  PLTSTRS. 


L 
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WRITEE 


CALL  WRITEE  (TIMEP) ; 

TIMEP  (*)  FL0AT  BIN.  A list  of  times  (seconds)  at  which 

temperature  distributions  are  to  be 
printed.  This  list  should  be  in 
ascending  order;  i.e.,  TIMEP(l)  < 
TIMEP(2)  < TIMEP(3)  etc. 

WRITEE  serves  the  same  function  as  WRITE  (Refs.  4 and  6), 
printing  out  temperature  distributions  (°F)  at  selected  times. 
WRITEE  suppresses  the  printing  of  fictitious  node  temperatures  if 
the  embedding-surface  node  numbers  read  by  READSU  are  positive  but 
prints  those  temperatures  if  the  embedding-surface  node  numbers 
read  by  READSU  are  negative.  After  printing  the  temperature  disti- 
butlon  of  the  body  being  analyzed,  WRITEE  prints  the  ablatlng-sur- 
face  temperatures  (°F)  and  recession  distances  (in.)  and  rates 
(in./s)  in  the  A direction  for  each  embedding  column,  distinguished 
in  the  output  by  its  embedding-surface  node  number. 

Each  temperature  distribution  listing  has  the  following 
general  appearance: 

25  - 1025.3  26  - 1040.9 

This  means  that  the  temperature  at  node  25  is  1025. 3°F  and  the 
temperature  at  node  26  is  1040. 9°F.  The  ablating  surface  listing 
has  the  appearance 

1 - 7055.2  4. 562E-04  3.411E-02 

for  each  embedding  column.  The  numbers  signify  that  the  ablating 
surface  in  the  embedding  column  with  embedding  surface  node  1 has 

a temperature  of  7055. 2°F,  is  located  4.562  x 10  ^ in.  from  the 
inertial  reference  surface  (see  Fig.  1)  in  the  A direction,  and  is 

moving  at  a rate  of  3.411  x 10  in./s  in  the  A direction  away 
from  the  inertial  reference  surface.  There  is  no  node  number  asso- 
ciated with  the  ablating  surface,  each  point  of  Intersection  be- 
tween the  ablating  surface  and  an  embedding  column  is  identified 
by  both  WRITEE  and  ABAER0E  by  the  associated  embedding-surface  node 
number . 
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When  WRITEE  prints  temperature  end  eblatlng-surface  infor- 
mation, it  prints  the  time  at  which  these  data  apply  before  print- 
ing anything  else.  After  everything  is  printed,  WRITEE  causes 
STEPE  to  print  the  time  step  used  and  the  maximum  allowable  time 
step,  based  on  the  temperatures  printed.  STEPE  also  prints  a line 
of  asterisks  to  separate  the  current  temperature  distribution  from 
the  next  one  to  be  printed. 

When  TIMS  is  greater  than  the  largest  value  of  TIMEP,  and 
there  is  more  than  one  value  of  TIMEP,  WRITEE  prints  temperature 
and  ablating  surface  Information  at  every  time  level.  If  TIMEP (1) 

< 0,  WRITEE  prints  information  at  every  time  step.  If  TIMEP (1) 

■ ±999,  additional  information  useful  for  debugging  is  printed. 

If  there  is  only  one  value  of  TIMEP,  TIMEP (1)  is  treated  as  a print- 
time Interval  if  it  is  positive,  and  WRITEE  prints  information  every 
TIMEP (1)  seconds.  If  there  is  more  than  one  value  of  TIMEP,  and 
TSTART  (saved  by  SET)  is  greater  than  the  first  several  values  of 
TIMEP,  no  printing  is  done  by  WRITEE  for  these  times. 

WRITEE  must  follow  ABAER0E  but  may  either  precede  or  follow 
PLTSTRS  (described  next). 
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PLTSTRS 


CALL  PLTSTRS  (PLTH,  #CAPS,  A,  B) 


(*)  FL0AT  BIN.  A list  of  times  (seconds)  at  which 
PLTSTRS  Is  to  save  data  for  plotting 
or  thermal-stress  purposes. 


PLTM 


FIXED  BIN.  The  number  of  nodes  In  the  flnlte-dif 
ference  model. 


#CAPS 


N.  A list  of  #CAPS  A coordinates. 

A(I)  - node  I's  A coordinate.  The  B or 
C coordinates  can  also  be  used  here. 


!N.  a list  of  #CAPS  B coordinates. 

B(I)  - node  I's  B coordinate.  The  A or 
C coordinates  can  also  be  used  here. 


PLTSTRS  Is  an  optional  procedure  that  can  be  used  for  sav- 
ing data  for  the  Rublnstein-Kleln  (Ref.  16)  plotter  code,  thermal 
stress  calculations,  and  restarting  the  thermal  analysis  computa- 
tions. PLTSTRS  has  six  files  (two  Input  and  four  output)  available 
for  these  purposes.  All  of  the  files  are  card  files  with  a logical 
record  length  of  80.  The  temperature  field  should  depend  on  only 
one  or  two  spatial  variables. 

JIMR0G  (named  after  APL's  J.  T.  Stadter  and  R.  0.  Weiss, 
whose  need  for  temperature  data  in  a special  format  for  thermal 
stress  computations  first  motivated  the  writing  of  PLTSTRS)  is  an 
Input  file  used  to  supply  comments  for  use  in  the  four  output 
files.  Up  to  twelve  lines  of  comments,  each  line  being  no  more 
than  forty  characters  long,  can  be  read  from  JIMR0G.  Only  the 
first  forty  columns  of  each  card  should  be  used.  There  has  to  be 
a DD  card  for  JIMR0G. 


NRPL0TS  is  an  output  file.  For  each  value  of  PLTM  and  all 
times  larger  than  the  maximum  PLTM,  PLTSTRS  stores  the  time  (sec- 
onds), the  A and  B coordinates,  temperature  (*F),  and  index  of 
each  real  node  and  uses  the  entire  contents  of  JIMR0G.  Fictitious 
node  data  are  also  stored  if  the  embedding  surface  node  nua&ers 
read  by  READSU  are  negative.  In  addition  to  node  data,  PLTSTRS 
saves  the  A and  B coordinates  and  temperatures  (*F)  of  the  ablating 
surface's  intersection  points  with  each  embedding  column  and  iden- 
tifies these  points  by  minus  the  embedding  surface  node  .lumbers  in 
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NRPL0TS . There  Is  one  card  Image  generated  for  each  node  In  a for- 
mat suitable  for  both  thermal  stress  calculations  and  use  with  the 
Rubinstein-Kleln  plotter.  If  there  Is  no  DD  card  for  NRPL0TS,  the 
file  is  Ignored. 

STRESS  Is  an  output  file  very  similar  to  NRPL0TS,  differing 
only  in  that  it  is  always  used  to  store  both  real  and  fictitious 
node  data  and  it  does  not  store  any  ablatlng-surface  temperatures. 
If  there  is  no  DD  card  for  STRESS,  the  file  is  ignored. 

RECPL0T  is  an  output  file  used  to  store  times  (seconds),  A- 
dlrection  recession  distances  (in.)  and  rates  (in./s),  and  embed- 
ding-surface node  numbers  for  each  embedding  column  at  times  cor- 
responding to  each  value  of  PLTM  and  greater  than  the  maximum  value 
of  PLTM.  These  data  can  be  used  with  the  Rubinstein-Kleln  plotter. 
RECPL0T  uses  the  first  line  in  the  JIMR0G  file  as  a descriptive 
title.  If  there  is  no  DD  card  for  RECPL0T,  the  file  is  ignored. 

NEXTRUN  is  an  output  file  that  saves  the  final  temperature 
distribution  and  time  for  use  in  a subsequent  job  or  job  step.  It 
uses  the  entire  contents  of  JIMR0G,  and  its  formatting  is  the  same 
as  that  used  for  NRPL0TS  and  STRESS.  If  there  is  no  DD  card  for 
NEXTRUN,  the  file  is  Ignored. 

PREVRUN  is  an  input  file  containing  the  contents  of  the 
NEXTRUN  file  of  a previous  job  or  job  step.  When  there  is  a DD 
card  for  PREVRUN,  PLTSTRS  overrides  SET  and  initializes  the  time 
and  temperatures  as  they  are  read  from  the  PREVRUN  file.  The  final 
time  in  this  case  is  taken  to  be  the  maximum  value  of  PLTM.  All 
this  is  done  on  the  first  pass  through  the  TACL.  If  there  is  no 
DD  card  for  PREVRUN,  the  file  and  the  actions  associated  with  it 
are  ignored. 

If  PLTM(l)  is  less  than  the  starting  time,  PLTSTRS  will  stop 
program  execution.  All  of  the  output  files  used  by  PLTSTRS  should 
be  allocated  to  sequential  data  sets  if  they  are  used  at  all. 
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8.  DEBUGGING  AIDS  AND  USEFUL  EXTERNAL 

IDENTIFIERS 

As  it  Is  with  other  SHTP  procedures,  the  SHTP-E  procedures 
hsve  several  external  identifiers  that  the  user  need  not  know. 
However,  there  are  some  external  identifiers  that  a user  may  de- 
sire to  use  in  his  main  program;  the  purpose  of  this  Section  is 
to  describe  them.  Each  identifier  will  be  introduced  by  its  DCL 
statement,  which  will  be  followed  by  a description. 


DCL  (DBGEB  INIT  ('0'B))  BIT(l)  EXT; 

DBGEB  is  a bit  flag  used  to  print  debugging  messages  throughout  the 
SHTP-E  procedures.  These  messages  are  printed  when  DBGEB  - 'l'B. 

A job  using  DBGEB  ■ 'l'B  should  be  run  for  only  very  few  time  steps 
because  of  the  copious  output  produced. 


DCL  (DBPTAB  INIT  ('O'B))  BIT(l)  EXT; 

When  DBPTAB  “ 'l'B,  ABAER0E  prints  the  constant-pressure  thermo- 
chemical  tables  that  it  generates  from  the  raw  data  obtained  from 
READABN.  This  is  done  only  when  ABAER0E  prints  ablating-surface 
information  as  dictated  by  its  print-time  argument. 


DCL  XCAP(*,7)  FL0AT  BIN  CTL  EXT; 

XCAP  is  the  array  of  geometric  parameters  read  by  READCPE.  XCAP(I,1) 
through  XCAF(I,7)  are  the  numerical  data  entries  read  in  columns  10 
through  79  on  the  cards  in  READCPE 's  geometry  file. 

DCL  T(*)  FL0AT  BIN  CTL  EXT; 

T(I)  is  the  temperature  of  node  I in  Rankine  degrees. 


DCL  Q(*)  FL0AT  BIN  CTL  EXT; 


Q(I)  is  the  sum  of  internal  source  heating  and  explicit  internodal 
heating  (in  Btu/h)  for  node  I.  For  each  time  level,  Q should  not 
be  altered  after  calling  CAPC0NE. 
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DCL  TIMS  FL0AT  BIN  EXT; 

TIMS  Is  the  current  thermal  model  time,  In  seconds. 


DCL  PP(*)  P0INTER  CTL  EXT; 


PP  Is  an  array  of  addresses  for  thermal  properties.  Its  dimension 
is  five  times  the  number  of  material  files  read  by  READPR,  when  PP 
is  used  to  access  a particular  property,  a based  variable  with  the 
declaration,  DCL  PRDUM(IOO)  BASED (PR)  FL0AT  BIN;  is  needed. 


,th 


For  the  J “ data  listing  of  material,  I,  PP,  and  PRDUM  can 
be  used  to  retrieve  the  properties  listed  below. 


PP(5*I-4)  - 


PP(5*I-4) 
PP(5*I-3) 
PP (5*1-2) 
PP (5*1-1) 
PP(5*I) 


PRDUM(l) 

Number  of  entries  in  property  list,  maxi- 
mum value  of  J. 

PRDUM (J+l) 

Temperature  (°R). 

PRDUM(J) 

Volumetric  specific  heat  (Btu/ft^-°R) . 

PRDUM(J) 

Thermal  conductivity  (Btu/h-f t-°R) . 

PRDUM(J) 

Emlssivity  or  gap  radiation  factor. 

PRDUM(J) 

Surface  absorptivity. 
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Computer  libraries  for  PL/ I source  code,  PL/ I optimizer-com- 
piler cross-referenced  listings,  and  compiled  and  link-edited  load 
modules  for  the  SHTP-E  procedures  have  been  created  under  the  Index 
BBE. RANDALL. 

BBE. RANDALL. EMBEDPL I contains  the  PL/ 1 source  code  for  all 
the  procedures  needed  to  support  SHTP-E.  Included  are  all  of  the 
procedures  described  In  Section  7.  It  also  contains  the  source 
code  for  DIKI,  BACKSUB,  and  SURFACE  which  are  called  by  ABAER0E; 

BAPSUB  which  Is  called  by  STEPS;  READRKD  which  Is  called  by  READPR; 
and  DECIDE  (Ref.  3)  and  PIF1D  (Ref.  3)  which  are  two  general-purpose 
Interpolation  subroutines.  The  library  member  names  for  the  source 
code  are  the  procedure  names;  e.g.,  ABAER0E's  source  code  Is  con- 
tained In  BBE. RANDALL. EMBEDPLI(ABAER0E). 

In  addition  to  the  procedure  source  codes,  BBE. RANDALL. EMBEDPLI 
has  three  other  members:  ENTRIES,  C0MMENTS,  and  LINKS. 

ENTRIES  contains  the  DCL  statements  for  all  of  the  procedures 
described  In  Section  7.  C0MMENTS  contains  a main  program  heading 
that  can  be  used  with  ENTRIES  to  avoid  blocksize  difficulties  when 
one  attempts  to  compile  a main  program  from  card  input.  Compiling 
a main  program  from  card  input  is  done  by  arranging  the  P.SYSIN 
cards  In  the  following  way: 

/ /P . SYSIN  DD  DSN-BBE . RANDALL . EMBEDPLI (C0MMENTS ) , DISP-SHR 

//  DD  * 

User  supplied  source  code. 

//  DD  DSN-BBE. RANDALL. EMBEDPLI (ENTRIES), DISP-SHR 

//  DD  * 

More  user  supplied  source  code. 

Use  of  C0MMENTS  AND  ENTRIES  In  this  way  makes  it  unnecessary 
for  the  user  to  write  DCL  statements  for  the  embedding  procedures 
he  has  to  call.  If  the  user's  source  code  is  in  a data  set  (e.g., 

BBE. USER. SOURCE)  with  a blocksize  of  6160  or  larger,  the  following 
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arrangement  can  be  used. 

//P.SYSLIB  DD  DSN-BBE. RANDALL. EMBEDPLI , DISP-SHR 
//P.SYSIN  DD  DSN-BBE. USER. S0URCE, DISP-SHR 
The  data  set  BBE. USER. SOURCE  should  contain  the  statement 
ZINCLUDE  SYSLIB (ENTRIES); 

In  the  calling  program. 

LINKS  is  a set  of  linkage  editor  Input  cards  that  can  be 
used  to  take  advantage  of  the  linkage  editor  overlay  feature  when 
an  executable  load  module  Is  created.  These  cards  can  be  accessed 
by  using  the  JCL  card 

//L.SYSIN  DD  DSN-BBE. RANDALL. EMBEDPLI (LINKS), DISP-SHR 

after  the  //L.SYSLM0D  statement  in  the  job  used  to  create  the  exe- 
cutable load  module. 

A listing  of  the  contents  of  any  or  all  of  the  members  of 
BBE. RANDALL. EMBEDPLI  can  be  obtained  by  using  the  LISTSRC  or  LISTPCH 
utilities  described  in  Ref.  17. 

Listings  of  the  PL/ I procedures  cross  referenced  by  the  PL/ I 
optimizer  in  BBE . RANDALL . EMBEDPLI  are  kept  in  members  of  a parti- 
tioned data  set  designated  by 

BBE . RANDALL . EBXREF (Procname) 

where  Procname  is  the  name  of  one  of  the  procedures.  The  CSDS 
utility  program  (Ref.  18)  can  be  used  to  obtain  these  listings. 

For  example,  to  obtain  a hard  copy  of  ABAER0E's  PL/ I cross  refer- 
enced listing,  one  can  use  the  following  JCL  cards. 

//LIST  EXEC  PGM-CSDS, TIME- (0,10) 

//SYSPRINT  DD  DUMMY 


Ref.  17.  "System  360  Model  91  User's  Guide,"  APL/JHU 
BCS-1:40,  Nov  1973. 

Ref.  18.  G.  L.  Bennett,  J.  C.  Hagan,  and  D.  C.  Tantlno, 
"Evaluating  Aeroshell  Materials  for  the  MJS/Multihundred  Watt  Heat 
Source,"  APL/JHU  ANSP-M-14,  Jul  1976. 
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//INPUT  DD  DSN-BBE . RANDALL . EBXREF (ABAER0E) .DISP-SHR 

u 

//0UTPUT  DD  SYS0UT-A,DCB*(RECF»^VBA,LRECL-137,BLKSIZE-1922) 

_ , j i 

A microfiche  copy  is  obtainable  by  replacing  SYS0UT-A  with  SYS0UT-Q 
on  the  //0UTPUT  card. 

Compiled  and  link  edited  SHTP-E  procedures  have  been  stored 
in  the  library  BBE. RANDALL. EMBEDLIB  whose  member  names  are: 

ABAER0E,  APFAC.  BACKSUB,  BAPSUB,  CAPC0NE,  DECIDE*,  DIKI,  PIF1D  , 

_ * * 

PLTSTRS , READ ABN,  READCPE,  READPR,  READRKD,  READSU,  READTM  , SET  , 

STEPE,  ST0RE,  SURFACE,  and  WRITES.  The  members  with  asterisks  were 
copied  directly  from  BBE. FRAZER. BASIC02.  BBE. RANDALL. EMBEDLIB  is 
needed  in  the  0L,  0PL,  or  0PLG  statement  when  an  executable  load 
module  is  created,  and  is  accessed  as  follows: 

//STEP  EXEC  0PL, LIB- 'BBE. RANDALL. EMBEDLIB ',.. . 
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10.  EXAMPLE  PROBLEM 


The  EM  procedures  have  been  used  successfully  in  reentry 
thermal  analyses  of  a General  Electric  Multihundred  Watt  Generator 
Heat  Source  Assembly  (MHW/HSA)  full-sphere  assembly  (FSA)  (Ref.  18) 
and  a Mound  Laboratory  radioisotope  heater  unit  (Ref.  19).  A 
specific  FSA  reentry  analysis  will  be  described  here  to  show  how 
the  EM  procedures  are  used. 

The  FSA  consists  of  a PuC>2  fuel  sphere  surrounded  by  an 

Iridium  post-impact  containment  shell  (P ICS),  that,  In  turn,  is 
surrounded  by  a T-50  wound  graphite  impact  shell  (G1S).  The  fuel 
sphere  Is  1.498  in.  in  diameter.  The  inside  diameter  of  the  PICS 
is  1.578  in.  and  it  has  a thickness  of  0.020  in.  The  Inside  diam- 
eter of  the  G1S  is  1.646  in.  and  its  thickness  is  0.460  in.  In  the 
example  given  here,  the  MHW/HSA  breaks  up  at  11.06  s into  a reentry 
flight  and  bore  FSA's  fall  from  this  point  to  Impact  at  191.66  s. 
Breakup  occurs  at  an  altitude  of  117  443  ft  at  a flight  path  angle 
of  -47.14°  at  a speed  of  29  830  ft/s.  The  GIS's  exterior  surface 
temperature  is  initially  2599°R  and  the  FSA's  temperature  distri- 
bution is  steady  at  breakup  time.  The  fuel-PICS  gap  contains  heli- 
um and  the  PICS-GIS  gap  is  evacuated. 

A cross  section  of  the  FSA  thermal  model  used  is  shown  in 
Fig.  4.  The  true  model  can  be  visualized  by  rotating  the  cross 
section  about  the  external  flow  field’s  axis  of  symmetry.  It  is 
assumed  that  the  temperature  distribution  is  axisymmetrlc  so  that 
it  depends  on  time,  radius  from  the  fuel  center,  and  angle  from 
the  forward  stagnation  point. 

The  FSA  analysis  was  carried  out  using  thermal  properties, 
thermochemical  data,  and  heating-rate  and  pressure-distribution 
data  obtained  from  various  sources.  PL/I  procedures,  JCL,  and  data 
needed  to  perform  the  analysis  are  stored  in  the  partitioned  data 
set  BBE. RANDALL. FSAEX,  a listing  of  which  is  obtainable  using  the 
LISTSRC  utility  (Ref.  16).  Table  5 lists  and  describes  the  members 
of  BBE. RANDALL. FSAEX. 


Ref.  19.  "Safety  Analysis  Report  for  the  MJS-Radioisotope 
Heater  Unit,"  Report  MLM-ML-76-48-0003,  Mound  Laboratory,  MlamlS' 
burg,  OH,  Aug  1976. 
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Table  5 

BBE.RANDALL.FSAEX 


Member 

Description 

GE0PLI 

PL/I  source  code  for  generating  geometric  data  for  the 
contents  of  the  members  GE0DATA,  SURFACE,  ABC00RDS, 

RSNPHIS,  and  0NEDS . 

GE0DATA 

Node  geometry  data  read  by  READCPE  for  the  thermal  model 
sketched  In  Fig.  4. 

SURFACE 

Embedding-surface  geometry  data  read  by  READSU. 

ABC00ROS 

Node  distance  from  reference-surface  data  read  by  READSU. 

RSNPHIS 

List  of  node  numbers  and  radial  and  angular  spherical 
coordinates  of  nodes  needed  by  PLTSTRS. 

0NEDS 

Node  geometry  data  for  a purely  radially  dependent  FSA 
temperature  distribution,  used  to  find  Initial  steady 
state. 

FUELTH 

Thermal  property  list  for  Pu02  fuel. 

GFPTH 

Thermal  property  list  for  Pu02-PICS  gap. 

GPSTH 

Thermal  property  list  for  PICS-GIS  gap. 

T50TH 

Thermal  property  list  for  GIS. 

IRIDIUM 

Thermal  property  list  for  PICS. 

ICPLI 

PL/I  source  code  for  finding  Initial  radially  dependent 

FSA  using  steady-state  procedure  described  In  Ref.  20. 

STJ0B 

JCL  necessary  to  execute  PL/I  procedure  In  member  ICPLI. 

IC2599 

Steady-state  temperature  distribution  generated  by 
steady-state  program. 

Ref.  20.  R.  K.  Frazer,  "An  Alternate  Solution  Technique  for 

Steady  State  Temperature  Analyals,"  APL/JHU  EM-4505,  Rev.  1,  26 

Oct  1973. 
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Table  5 (cont'd) 


Member 

Description 

FC0NV 

List  of  convective  heating  rate  factors  used  to  define 
convective  heating  distribution  on  the  FSA's  outer  sur- 
face. 

TRAJ 

Trajectory  data  generated  by  program  described  in  Ref. 

15. 

ANPLI 

PL/ I source  code  for  the  main  program  using  the  SHTP-E 
procedures  to  solve  the  FSA  reentry  problem. 

LIBCR 

JCL  for  creating  a compiled  and  link-edited  library 
using  the  contents  of  member  ANPLI. 

LIBREP 

JCL  for  replacing  the  libary  created  by  the  job  in  mem- 
ber LIBCR  with  a revised  library. 

CHANGE 

JCL  for  creating  an  executable  load  module  using  the 
contents  of  BBE . RANDALL . EMBEDLIB  and  the  library  created 
by  the  job  in  member  LIBCR  or  revised  by  the  job  in  mem- 
ber LIBREP. 

J0B1 

JCL  for  a job  using  the  executable  load  module  created 
by  the  job  in  member  CHANGE  in  the  FSA  thermal  analysis 
beginning  with  the  initial  FSA  steady-state  data  in  mem- 
ber IC2599.  The  JCL  shows  how  to  create  output  data 
sets  from  the  contents  of  the  output  files  used  by  the 
procedures  ABAER0E  and  PLTSTRS. 

J0BM1 

JCL  for  a job  that  begins  with  the  data  created  at  mem- 
ber J0Bl's  last  time  level.  The  JCL  shows  how  to  revise 
the  data  sets  created  by  the  job  in  member  J0B1. 

The  aerothermochemical  data  were  taken  from 


BBE. FRAZER. MHWDATA(TCH14X39) . 

A PL/I  optimizer-compiler  cross-referenced  listing  of  the 
contents  of  BBE. RANDALL. FSAEX(ANPLI)  is  kept  in  the  sequential  data 
set  BBE . RaNDALL . FSAXREF . The  printed  outputs  from  the  jobs  listed 
in  BBE. RANDALL. FSAEX(J0B1)  and  BBE. RANDALL. FSAEX(J0BM1)  are  kept  in 
the  sequential  data  sets  BBE . RANDALL . FSAPR1  and  BBE . RANDALL . FSAPRM1 , 
respectively.  Listings  of  the  contents  of  these  three  sequential 
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data  sets  can  be  obtained  using  the  JCL  given  In  Section  9 for  list- 
ing the  contents  of  BBE . RANDALL. EBXREF(ABAER0E) . 

A microfiche  or  hard-copy  listing  that  Includes  PL/ I cross- 
referenced  listings  of  all  the  embedding  procedures,  the  members 
of  BBE. RANDALL. FSAEX  discussed  above,  the  aerothermochemlcal  data, 
a PL/I  cross-referenced  listing  of  the  main  FSA  thermal-analysis 
procedure,  and  the  two  printed  outputs  kept  In  BBE . RANDALL . FSAPR1 
and  BBE. RANDALL. FSAPRM1  can  be  obtained  using  the  following  JCL 
and  data  cards. 

//Job  statement 

//LIST  EXEC  PGM-CSDS, TIME- (0,10) 

//SYSPRINT  DD  DUMMY 
//INPUT  DD  DATA 

//J0BNAME  J0B  (ACCT//,C,U,N) , "Title" 

//  'USER  ID',RD-R 
//*R0UTE  Destination 
/* 

//  DD  DSN-BBE. RANDALL. FSAEX (Member ),DISP-SHR 
//0UTPUT  DD  SYS0UT-J, 

//  DCB- (RECFM-FB .LRECL-80 , BLKSIZE-6160) 

/* 

The  above  JCL  contains  a Job  that  submits  another  job.  The 
file  INPUT  is  a concatenation  of  a user-supplied  job  statement  and 
destination  card  with  JCL  necessary  to  produce  the  desired  output. 
Using  SYS0UT  - J,  the  file  0UTPUT  submits  the  Job  assembled  in 
INPUT.  The  JCL  of  this  Job  is  printed  at  remote  terminal  destina- 
tion (e.g.  RM001),  but  the  printed  (Member  - PRXREF)  or  microfiche 
(Member  - FICHXREF)  output  is  produced  at  the  APL  computing  center 
(Building  3). 
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NOMENCLATURE 


English 
A,  B,  C 


Orthogonal  coordinates,  m ox  In. 

2 2 

Embedding  surface  node  area,  m or  ft 


Lagrange  Interpolation  coefficient  for  temperatures 


Heat  flux  weighting  factor 


b(Ei,Ej) 


Lagrange  Interpolation  coefficient  for  finding 
ablating  surface  gradient,  m ^ or  ft  ^ 

Logorlthmic  ablation- to-dlf fusion  mass-loss  ratio 

Thermal  capacitance,  J/K  or  Btu/°R 


Heat  transfer  Stanton  number 


Mass  transfer  Stanton  number 


Lagrange  Interpolation  coefficient  for  finding 
ablating  surface  temperature  gradient,  m or  ft 

Distance  of  node  1 In  embedding  column  from  Iner- 
tial reference  surface,  m or  in. 

Distance  of  embedding  surface  from  fixed  reference 
surface,  m or  In. 

Recession  distance  iteration  parameter 


H<T.) 


Functional  relationship  describing  ablating  sur- 
face temperature  dependence  of  the  ratio  of  me- 
chanical erosion  mass-loss  rate  to  thermochemical 
mass-loss  rate 


Thermal  conductance  for  heat  flow  between  nodes  1 
and  J,  W/K  or  Btu/h-#R 


- 87  - 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

LAUREL  MARYLAND 


k 

* 

Le 

m 

P 

QF 

Q 

qchem 

qcond 

q 

nconv 

qEXT 

qi 

qradin 

qradout 

qrem 


Thermal  conductivity  of  ablating  material,  W/m-K 
or  Btu/h-f t-*R 

Equivalent  conduction  path,  m or  ft 
Lewis  number 

2 

Mass  loss  flux  due  to  ablation,  kg/m  -8  or 
lbm/ft^-s 

Tabular  pressure  In  thermochemical  table;  local 
pressure,  atm 

Fictitious  heating  rate  applied  to  embedding  sur- 
face, W or  Btu/s 

Sum  of  distributed  heat  source  and  explicit  inter- 
nodal  heating  rates,  W or  Btu/h 

Ablating  surface  heat  flux  due  to  chemical  reaction, 
W/m^  or  Btu/ft^-s 

Heat  flux  conducted  frr.m  ablating  material  to 

2 2 

ablating  surface,  W/m  or  Btu/ft  -s 

External  convective  heat  flux  to  ablating  surface, 
W/m^  or  Btu/ft^-s 

Total  external  heat  flux  to  ablating  surface, 

W/m^  or  Btu/ft^-a 

Internal  heat  generation  rate  for  node  1,  W or 
Btu/h 

Incident  radiative  heat  flux  on  ablating  surface, 
W/m^  or  Btu/ft^-s 

Radiation  relief  from  ablating  surface  to  space, 
W/m^  or  Btu/ft^-s 

Energy  flux  required  to  remove  ablating  material, 
W/m^  or  Btu/ft^-s 


0 


( 


J 


.1 


.1 


j 


] 


I 


t 

I 

I 

I 

f 


i j 
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Greek 

e 


AT 


0 

L 


At 


Sum  of  explicit  conductances  for  node  1,  W/K  or 
Btu/h-#R 

Sum  of  Implicit  conductances  for  node  1,  W/K  or 
Btu/h-°R 

A-dlrection  recession  distance,  m or  in. 

Temperature,  K or  °R 

Adiabatic  thermal  response,  K or  °R 

Unit  incremental  thermal  response,  K/W  or  °R-s/Btu 

Temperature  used  to  estimate  ablating  surface 
temperature  gradient,  K or  °R 

Time,  s or  h 

B-directlon  temperature  approximation,  K or  °R 
C-direction  temperature  approximation,  K or  °R 
Surface  thermochemical  table  dependent  variable 


Ratio  of  actual  ablation  mass-loss  rate  to 
diffusion  limited  mass-loss  rate,  uncorrected  for 
blowing 

Blowing-corrected  ratio  of  actual  ablation  mass- 
loss  rate  to  diffusion-limited  mass-loss  rate 

Kronecker  delta:  * 1,  if  i ■ E;  6iE  = 0,  if 

i + E 

Ablating  surface  temperature  change  during  material 
re-oval;  temperature  change  during  rough  search  for 
ablating  surface  temperature,  K or  °R 

Time  step,  t^  - t^,  s or  h 
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e’(T) 


Subscripts 


I.  J,  K 


2 2 

Newton  Raphson  error  function,  W/m  or  Btu/ft  -s 
for  embedded  ablating  surface,  K or  °R  when  ablat- 
ing and  embedding  surfaces  coincide 

de/dT 

Blowing  parameter 

Cosine  of  angle  between  A direction  and  ablating 
surface  normal 

Lagrange  interpolation  polynomial  approximation 

of  3T,V3A,  K/m  or  °R/ft 
A 

Lagrange  interpolation  polynomial  approximation 
of  STy/SA,  K/W-m  or  s-°R/Btu-ft 

Logarithmic  tabular  pressure  in  thermochemical 
tables;  local  logarithmic  pressure 

3 3 

Density  of  ablating  material,  kg/m  or  lbm/ft 
Sum  of  temperature  coefficients,  J/K  or  Btu/°R 


Blowing  corrected  value 

Embedding  surface  node 

Gas  flowing  over  ablating  surface 

Three  successive  nodes  in  embedding  column 

Node  being  analyzed;  tabular  pressure  entry  in 
thermochemical  tables 

Neighboring  node;  tabular  mass-loss  ratio  entry 
in  thermochemical  tables 

Mechanical  erosion 

Ablating  surface 

Ablating  surface  intersection  with  embedding 
column  having  embedding  surface  node  E 
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TC 

Thermochemical  ablation 

Superscripts 

A,  B,  C 

A-,  B-,  C-direction  implicit  internodal  links 

E 

Explicit  internodal  link 

l 

Time  level  index,  time  step 

+ 

Time  level  4+1.  No  superscript:  time  level 

n 

Iteration  number 

Operators 

Difference  operator,  AijTi  “ Ti  “ Tj 

VBC 

Gradient  in  B-C  plane,  m ^ or  ft  ^ 

n 

Product 

z 

Summation 

Miscellaneous 

F 

Vector  of  all  values  of  F for  all  nodes 

Acronyms 

AF 

Approximate  factorization 

EBE 

Energy-balance  equation 

EM 

Embedding  method 

FD 

Finite  difference 

SEBE 

Ablatlng-surface  energy-balance  equation 

SHTP 

Standard  Heat  Transfer  Program 

SHTP-E 

Embedding  Standard  Heat  Transfer  Program 
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STT 

TACL 


Surface  thermocheaical  table 
Thermal-analysis  computation  loop 


